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CFD  OF  COMPLEX  THREE-DIMENSIONAL 
MULTIPHASE  FLOWFIELDS 

_ PHASE  HSBIR  Topic  No.  A95-103 


1.0  INTRODUCTION 

The  post-hit  fate  of  injurious  chemical  or  biological  (C/B)  threat  vehicle  payloads  is  a 
critical  element  of  missile  lethality  programs.  Substantive  testing  and  varied  modeling  efforts  have 
been  performed  in  support  of  missile  C/B  lethality.  The  physics  and  thermochemical/rheological 
aspects  of  the  post-hit  scenario  are  extremely  complex  and  hence,  substantive  uncertainties  remain 
in  estimating  lethality.  A  scenario  of  primary  concern  has  been  that  of  higher  altitude  missile 
intercept.  At  higher  altitudes,  bulk  liquid  from  ruptured  payloads  is  initially  exposed  to  low 
pressures  and  mechanisms  of  cavitation  and  vaporization  can  lead  to  its  partial  disintegration  (with 
gas/bubble  loading  and  volatility  characteristics  serving  as  key  parameters).  As  the  remaining  bulk 
or  droplets  head  towards  earth,  aerodynamic  breakup  in  entering  the  atmosphere  has  been  found  to 
play  a  discriminating  role,  with  the  key  features  schematized  in  Figure  1.1.  The  fate  of  surviving 
droplets  has  also  been  found  to  have  a  strong  dependence  on  atmospheric  properties  (e.g.  humidity) 
with  droplet  surface  microphysics  (e.g.  condensation)  found  to  effect  their  survivability. 

The  higher  altitude  engagement  scenario  requires  that  maximum  destruction  of  the  C/B 
payload  be  done  at  intercept.  Due  to  the  lack  of  sufficient  oxygen,  conventional  high  energy 
(chemical)  mechanisms  are  not  applicable  for  providing  a  more  robust  kill  and  for  post-hit 
neutralization  of  the  undamaged  payload.  Post-hit  C/B  agent  fate  remains  largely  dependent  on 
natural  processes  (cavitation/vaporization,  aerobreakup,  etc...)  with  little  available  in  the  form  of 
countermeasures  during  bulk  cloud  flyout  and  its  progression  through  the  atmosphere  to  the 
ground. 

Lower  altitude  engagement  scenarios,  such  as  cruise  missile  (or  endo-engagement  ballistic 
missile)  intercept,  permit  supplementing  kinetic  energy  with  chemical  (explosive)  energy,  both  on 
impact,  and  as  an  early-time  countermeasure.  Such  scenarios  have  not  been  investigated  in  any 
detail,  to  date.  What  has  been  investigated  in  some  detail  is  the  neutralization  of  C/B  agents  stored 
in  bunkers,  in  programs  such  as  Agent  Defeat.  Common  concerns  in  all  C/B  related  scenarios  are 
those  of  initial  (post-hit)  agent  characteristics  (whether  escaping  from  a  missile  or  a  bunker 
canister),  the  response  of  agents  to  their  often  highly  dynamic  surroundings,  and,  the 
thermochemical/rheological  properties  of  agents. 

Data  is  evidently  required  to  support  model  calibration/validation  and  the  models  are 
expected  to  provide  some  insight  into  what  is  occurring  in  varied  C/B  scenarios,  with  different 
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expectations  from  different  categories  of  models.  It  is  important  to  clarify  the  relationship  between 
models,  enabling  data,  and  predictive  expectations,  since  this  has  led  to  some  confusion  in  spanning 
the  chasm  between  research  developments  and  systems  requirements.  Table  1.1  provides  this 
categorization. 


Table  1.1.  C/B  Model  Categorization. 


MODEL  TYPE 

ENABLING  DATA 

ROLE/PREDICTIVE 

EXPECTATIONS 

Research  Model 

•  First-Principles  for 
resolvable  scales 

•  Generalized  capabilities  for 
modest  scenario 
complexity 

•  Unit  problem  laboratory  data 
with  detailed  diagnostics 

•  Property  data 

•  Micro-physics  data 

•  Shed  insight  into  physics 

•  Help  interpret  ground/flight 
test  data 

•  Guide  construction  of 
engineering  model 

Engineering  Model 

•  Scenario  specific 

•  Difficult  physics  is  tuned 

•  Ground/flight  test  data 

•  Research  model  solutions  of 
unit  problems 

•  Predictive  capabilities  for  range 
of  calibrated  scenarios 

•  Provide  sensitivities  to 
variations  from  calibrated 
scenarios 

Systems  Model 

•  End-to-end  scenarios 

•  Tuned  physics  or  lumped 
parameter 

•  Engineering  model  solutions 
of  scenarios 

•  Scenario  data  sets 

•  End-to-end  scenarios  or 
portions 

•  Bottom  line  lethality 
predictions 

Referring  to  Table  1.1,  a  key  point  is  that  there  should  be  a  hierarchy  of  model  development, 
with  research  models  supporting  construction  of  engineering  models,  and  engineering  models 
providing  solution  data  sets  for  systems  model  calibration.  Research  models  require  detailed  sets 
of  unit  problem  data  (e.g.  blob  breakup  data  with  high  resolution  imagery  of  the  surface  details  and 
of  the  stripped  droplet  ligaments),  and  should  be  initiated  using  fluids  for  which  properties 
(viscosity,  surface  tension,...)  are  well  known.  Microphysics  data  is  required,  for  example,  to 
support  simulating  a  spores  response  to  its  environment,  that  would  contain  detailed  property 
variations  within  the  spore.  Measurements  of  spore  survivability  percentages  after  exposure  to  high 
temperatures  for  varied  time  intervals,  in  a  controlled  environment,  does  not  provide  what  is  needed 
by  a  research  model.  Research  models  can  directly  support  ground/flight  tests  by  supplementing 
test  data  with  key  missing  parameters.  For  example,  in  missile/canister  ground  tests  where  post-hit 
droplet  cloud  data  has  been  obtained  (imagery,  IR...),  research  models  with  advanced  droplet 
methodology  can  be  utilized  to  estimate  droplet  characteristics  (size  on  expulsion,  etc.). 
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In  the  Phase  II  SBIR  effort  to  be  described,  research  multi-phase  Navier-Stokes  (NS) 
models  were  extended/formulated  to  provide  the  framework  for  analyzing  aerodynamic  breakup  in 
both  low  and  high  speed  slip  velocity  (V  [droplet]  -  V  [gas])  regimes  with  an  emphasis  on  bulk 
liquid  payloads.  The  initiating  process  is  the  deformation  of  the  bulk  liquid.  Our  structured  grid 
research  (NS)  model,  CRAFT,  also  permits  performing  neutralization  studies  as  it  contains 
generalized  gas-phase,  liquid  (droplet)  and  solid  (particulate)  combustion  capabilities  stemming 
from  earlier  research  efforts  related  to  gun,  rocket  and  airbreathing  propulsion.  Section  2 
summarizes  the  CRAFT  code  model  capabilities  relevant  to  bulk  liquid  post-hit  scenarios. 

A  detailed  flyout  and  aerobreakup  study  of  a  high  speed  blob  is  presented  in  Section  3, 
showing  the  complex  surface  deformation  characteristics  and  the  droplet  wake  cloud  structure 
behind  the  blob.  The  full  dynamics  from  blob  flyout  to  droplet  wake  cloud  formation  is  shown. 

In  Section  4,  we  examine  the  neutralization  of  droplet  clouds  by  conventional  chemical 
explosive  mechanisms.  The  studies  exhibit  the  marked  sensitivities  to  the  droplet  sizes  and  to  the 
timing  of  the  explosion  relative  to  the  clouds  trajectory.  With  sufficient  energy  and  correct 
timing/charge  locations,  the  possibility  of  neutralization  countermeasures  for  C/B  threats  for  endo- 
intercepts  can  be  made  viable. 

Section  5  describes  the  engineering  model,  SDROP,  which  tracks  a  droplet,  from  high 
altitude  to  the  ground,  accounting  for  trajectory,  deformation,  variable  temperature,  and  vaporization 
effects.  The  importance  of  including  these  effects  and  their  impact  on  droplet  survivability  are 
discussed.  Use  of  research  models  to  study  shape  deformation  of  droplets  is  also  addressed. 

In  the  course  of  performing  the  work  described,  limitations  in  the  ability  of  the  CRAFT  code 
to  analyze  varied  aspects  of  the  bulk  liquid  flyout/breakup  problem  with  high  accuracy  were 
identified  which  would  require: 

1  use  of  unstructured  grid  numerics  with  dynamic  adaption  features  to  highly  resolve  the 
gas/liquid  interface; 

2  formulation  of  multi-phase  numerics  which  can  concurrently  analyze  incompressible, 
viscous  (Newtonian  or  non-Newtonian)  liquids  and  highly  compressible  gases  over  a  broad 
range  of  speeds;  and, 

3  cavitation  capabilities. 

No  present  code  has  such  capabilities.  In  chapter  6,  we  describe  significant  progress  toward 
addressing  these  requirements  in  an  all-speed,  all-fluid  variant  of  the  CRUNCH  unstructured  grid 
NS  code. 
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DYNAMIC  STRIPPING  OF  DENSE  DROPLET  DILUTE  DROPLET  VAPOR  TRAIL 

INTERACTION  OF  DROPLETS/  CLOUD-  CLOUD  IN  FAR  RISING  DUE  TO 

BULK  LIQUID  WITH  FILAMENTS  SECONDARY  WAKE -FURTHER  BUOYANCY 

ATMOSPHERE-  FROM  LIQUID  BREAKUP  AND  VAPORIZATION  EFFECTS 

DEFORMATION.  OCCURS  AT  VAPORIZATION  IN 

ELONGATION  OF  GAS/LIQUID  HOT  WAKE 

FILAMENTS  INTERFACE 

Fig.  1.1.  Schematic  of  the  breakup  of  a  mass  of  bulk  liquid  at  high  relative  velocity. 
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2.0  CRAFT  GAS/LIQUID  CODE  OVERVIEW 


2.1  INTRODUCTION 

The  gas/liquid  multi-phase  formulation  in  the  CRAFT  code  treats  two  different  classes  of 
multi-phase  problems.  The  first  category  deals  with  heterogeneous  gas/liquid  mixtures  where  the 
liquid  is  a  bulk  continuum  fluid.  Problems  in  this  class  include  gas/liquid  jets  and  bulk  liquid 
flyout.  For  this  class  of  problems,  the  formulation  we  discuss  here  entails  the  solution  of  a  unified 
set  of  conservation  equations  for  the  generalized  gas/liquid  mixture.  The  interface  between  the  gas 
and  liquid  phases  is  captured  as  part  of  the  solution  procedure,  permitting  its  shape  to  evolve  in  a 
transient  fashion  as  it  interacts  with  the  gas.  The  conservation  equations  are  cast  in  strong 
conservation  form  using  a  compressible,  density  based  formulation.  An  upwind  flux  difference 
scheme  is  utilized  to  integrate  the  equations  which  captures  shock  waves  as  well  as  contact  and 
gas/liquid  interfaces.  Pressure  waves  can  propagate  across  inter-phase  boundaries  in  this  scheme 
with  minimal  numerical  diffusion,  and  with  speeds  consistent  with  the  thermodynamics  of  the  multi¬ 
phase  composition.  Both  the  gas  and  liquid  components  can  be  represented  using  generalized 
compressible  equations  of  state. 

The  second  class  of  is  that  of  a  generalized  continuum  fluid  which  contains  a  dense 
dispersed  phase  (particles  or  droplets)  which  is  in  non-equilibrium.  Problems  in  this  category 
encompass  a  wide  range  of  applications  but  for  C/B  problems,  they  encompass  the  analysis  of  the 
fate  of  droplet  or  particulate  clouds.  This  system  of  equations  retains  its  hyperbolic  characteristics 
with  a  dispersive  wave  system  which  is  consistent  with  transient  multi-phase  physics.  Furthermore, 
the  equation  system  contains  source  terms  proportional  to  the  porosity  gradient. 

2.2  EQUATION  SYSTEM  FOR  GENERALIZED  GAS/LIQUID  MIXTURE 
2.2.1  Gas/  Bulk  Liquid  Mixtures 

We  begin  by  describing  the  numerical  formulation  for  a  generalized  gas/liquid  mixture  in 
the  CRAFT  code.  Since  the  bulk  liquid  represents  the  source  of  mass  for  the  subsequent  post¬ 
breakup  dispersed  phase,  its  mass  may  dominate  the  flowfield  in  specific  regions  of  the  domain, 
and  computing  the  bulk  liquid  shape  as  accurately  as  possible  is  a  priority.  Full  simulations  of 
multi-fluid  problems  are  generally  performed  either  by  explicitly  tracking  the  position  of  the  gas 
liquid  interface  using  Lagrangian  markers  [1]  or  capturing  it  using  a  Volume  Of  Fluid  method 
(VOF)  [2].  While  these  methods  make  it  possible  to  ascertain  the  interface  location  explicitly,  they 
require  a  predetermined  logic  to  manage  cases  where  the  interface  is  very  convoluted.  The  third 
alternative  which  we  employ  is  that  of  interface  capturing  via  writing  the  NS  equations  for  a  mixture 
of  gas  and  liquid  in  strong  conservation  form  here,  the  interface  between  the  gas  and  the  liquid 
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phase  is  the  localized  volume  which  contains  cells  with  both  gas  and  liquid.  Via  use  of  a  low 
numerical  diffusion  upwind  scheme  to  integrate  the  NS  equations,  the  numerical  thickness  of  this 
interface  can  be  restricted  to  a  few  computational  grid  cells.  This  is  the  method  implemented  in 
CRAFT. 

The  gas/liquid  conservation  equations  are  given  as: 


dQ  dE  dF  dG 


(2.1) 


with  the  vector  of  independent  variables  Q  and  convective  flux  E  being  defined  as  follows: 


Q  = 


Pm 

pnum 

r  fn 

{pmumum+exp) 

PmVm 

(pmEmvm+iyP) 

PmWm 

;  E  = 

(pA^+tP) 

e* 

(em  +  P)Um 

Pi 

"  £ 

.  Pn~\  . 

P  ,t/ 

_  rn-l  m 

(2.2) 


The  right-hand-side  of  Eqn.  (2.1)  includes  the  effects  of  the  viscous  forces,  and  the  source 
term,  S,  contains  terms  related  to  phase  change  (due  to  breakup,  vaporization  or  burning)  and 
interactions  between  the  gas  phase  and  dispersed  phase.  The  viscous  terms  are  similar  to  those  for 
a  single-phase  simulation.  The  different  components  of  S,  which  represent  the  inter-phase  coupling 
terms  between  the  flow  solvers,  are  presented  in  detail  in  the  subsections  below. 

The  first  five  equations  solve  for  the  conservation  of  mixture  mass,  momentum,  and  energy. 
The  subsequent  equations  solve  for  the  (n-1)  gaseous  species  (the  gas  can  be  a  generalized  mixture 
and  the  equations  permit  gas  phase  combustion  and  inclusion  of  generalized  combustion  products). 
The  liquid  mass  is  obtained  from  the  global  mixture  continuity  using  the  following  relationship: 

Pl<Pl  =  Pm~  Pg<t>g  (2-3> 

Here,  <pL,  denotes  the  liquid  volume  fraction  while,  (pg ,  represents  the  gas  volume  fraction.  We  note 
that  at  a  given  grid  point,  a  single  velocity  is  utilized.  Since  the  gas/liquid  interface  is  well  resolved, 
this  does  not  pose  a  problem  since  it  enforces  local  phase  equilibrium  only  at  the  interface  region 
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where  gas/liquid  coexist  in  a  computational  cell.  The  more  important  issue  is  being  able  to  capture 
the  interface  in  as  few  cells  as  possible:  a  problem  akin  to  capturing  a  strong  shock  in  gas 
dynamics. 

The  equation  of  state  for  gas  species  can  be  that  of  a  perfect  gas  mixture,  or  for  high 
pressures,  can  be  specified  using  non-ideal  Virial  or  Noble- Abel  equations  of  state.  For  the  Virial 
equation  of  state,  the  partial  pressure  for  each  gaseous  species  is  defined  as 


p  -  MI 

M; 


[l  +  B„p  +  Cvp2] 


while  it  is  given  as 


9  .  MI 

*  Mi[l -bp] 


(2.4a) 


(2.4b) 


for  the  Noble- Able  equation  of  state. 

The  liquid  equation  of  state  can  be  specified  either  in  tabular  form  or  as  an  analytical 
function,  if  available.  For  simple  compressible  liquids,  we  generally  use  the  following  analytical  ex¬ 
pression 


P 

L  *2 


f  \k2 

<Po 


-l 


(2.5) 


where  p0  is  the  liquid  density  at  a  reference  state,  is  the  bulk  modulus  at  zero  pressure,  and  the 
index  k2  is  a  measure  of  the  liquid's  compressibility. 

The  enthalpy  for  the  generalized  gas/liquid  mixture  is  obtained  using  the  following 
thermodynamic  relation: 


dh  =  C (T,P)rdT+ 


/p  "  ar 


dP 


which  reduces  to  the  following  expressions  for  varied  limiting  situations: 

For  an  ideal  gas:  h  =  J Cp(T)dT 

T 

For  the  Noble- Abel  EOS:  h  =  J  CpdT + bP 

T 


(2.6) 


(2.7a) 


(2.7b) 
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For  a  compressible  liquid: 


h  =  \CJT)dT  +  - 


k-y-l 


+ 


P(* 2-1) 


(2.7c) 


To  compute  the  thermodynamic  variables  we  need  to  decode  the  temperature  and  volume 
fraction  of  liquid  from  the  conserved  variables,  Q.  This  is  done  by  iterating  on  the  following  two 
conditions  simultaneously: 

hm  -  —  =  em=  —  - 1  /  2(u2  j  (< condition  on  temp.)  (2.8a) 

Pm  Pm 


Pg-PL  ( condition  on  (j)g )  (2.8b) 

For  temperature  nonequilibrium  between  gas  and  liquid,  the  liquid  temperature  is  specified 
independently  and  the  equation  for  internal  energy  provides  us  with  the  gas  temperature. 

2.2.2  Gas/Bulk  Liquid  Mixture  with  a  Dense  Dispersed  Phase 

In  this  section,  we  briefly  present  extensions  to  the  generalized  gas/bulk  liquid  formulation 
given  above  to  account  for  a  dense  dispersed  phase  of  droplets  or  particulates.  The  droplets  are  in 
non-equilibrium  with  the  continuum  phase  and  can  occupy  a  significant  fraction  of  the  grid  cell 
volume.  The  vector  of  dependent  variables  Q  and  the  convective  flux  terms  in  Eqn.  (2.2)  for  the 
continuum  phase  gets  modified  as  follows: 
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Here,  the  term  am  denotes  the  fraction  of  the  total  cell  volume  occupied  by  the  continuum  phase 
and  is  defined  as 


a, 
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(2.10) 
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where  Vp  is  the  total  volume  occupied  by  the  dispersed  phase  in  the  cell,  which  is  required  to  be 
incompressible.  Furthermore,  the  continuum  phase  equations  have  additional  source  terms  arising 
from  volumetric  and  non-equilibrium  drag  effects.  The  source  term,  D,  is  written  as 
D  =  S  +  Fd  +  Fp,  and  includes  the  effect  of  three  independent  physical  phenomena.  The  source 
term,  S,  contains  the  mass  and  energy  transfer  terms  due  to  vaporization/combustion  of  droplets 
(and  to  gas  phase  kinetics  as  well).  The  source  term,  FD,  contain  the  inter-phase  drag,  while  FP 
contains  source  terms  due  to  volumetric  effects.  These  source  vectors  are  given  as  follows: 
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(2.11a) 


(2.11b) 


The  details  of  the  drag  coefficients  Ap  and  Bp  are  available  from  correlations  for  packed  beds  or 
from  correlations  for  rocket  propulsion  if  sufficiently  dilute  [3]. 

The  system  described  in  Eqns.  (2.1)  and  (2.9)  has  a  real  set  of  eigenvalues  associated  with 
it.  We  begin  by  addressing  the  numerical  characteristics  of  this  system  of  equations.  The 
eigenvalues  of  the  inviscid  flux  Jacobian  in  Eqn.  2.9  are  hyperbolic  and  have  the  acoustic  speed 
associated  with  the  carrier  phase.  This  is  the  so-called  "frozen"  acoustic  wave  which  is  initially 
transmitted  before  the  particles  have  reacted  to  the  change  in  the  carrier  fluid  conditions.  The  drag 
term  (source  term  Fd)  slows  the  continuum  velocity  down  and  eventually,  at  large  times,  we  obtain 
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the  "equilibrium"  acoustic  wave  which  is  associated  with  the  equilibrated  conditions.  The  source 
term  FP  provides  local  acceleration/retardation  to  the  gas  based  on  the  gradient  of  porosity. 

2.2.3  Lagrangian  Dispersed  Phase  Formulation 

The  equations  for  the  dispersed  phase  are  presented  here  in  a  Lagrangian  manner  for  mass, 
momentum,  and  energy  of  individual  numerical  droplets  or  particles.  The  porosity  of  the  dispersed 
phase  is  computed  from  summations  at  the  individual  locations  of  the  particles.  We  note  that 
computing  this  Eulerian  quantity,  without  spurious  oscillations  resulting  from  grid  dependency,  is 
critical  since  it  impacts  the  numerical  stability  of  the  scheme.  The  droplet/particle  equations  are 
given  as: 


dX 


2-=V 
dt  p 
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where  the  drag  and  heat  transfer  coefficients  are 
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(2.12a) 

(2.12b) 

(2.12c) 

(2.13) 

(2.14) 


where  ps  is  the  material  density  of  the  particulate  material,  Cd  is  the  drag  coefficient,  Cp  is  the 
specific  heat  of  the  gas,  and  Nu  the  Nusselt  number.  In  the  results  discussed  below,  Cd  and  Nu  are 
computed  using  correlations  derived  for  rocket  propulsion  applications  which  are  valid  over  a  wide 
range  of  lag  Reynolds  and  Mach  numbers  (see  Ref.  3). 

An  additional  equation  is  implemented  to  take  into  account  the  change  of  particle  radius  rp  in 
time,  as  a  result  of  evaporation  or  burning,  or  of  secondary  breakup. 


*jl 

dt 


=  ra+rb 


(2.15) 
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We  note  that  the  momentum  equation  as  shown  retains  an  intergranular  stress  term 
formulated  primarily  for  a  dense  solid-particulate  phase  which  represents  the  additional  stress 
experienced  by  the  solid-phase  below  the  settling  porosity  where  the  particles  touch  each  other. 
This  effect  may  not  be  appropriate  for  liquid  droplets  whose  phenomenology  is  different.  For  very 
high  loading,  liquid  droplets  will  coalesce,  and  become  a  continuum  fluid  whose  pressure  may  then 
be  described  by  an  equation  of  state  (as  in  Eqn.  2.5).  For  the  results  presented  here,  coalescence 
has  not  been  modeled  and  represents  an  area  of  future  work. 

2.2.4  Eulerian  Dispersed  Phase  Formulation 

A  new,  second-order  upwind  Eulerian  solver  was  incorporated  in  the  CRAFT  code  and  coupled  to 
the  gas/liquid  mixture  through  source  terms.  Conservation  equations  for  the  Eulerian  continuum 
cloud  (written  in  dilute  form  with  no  volumetric  effects)  are  solved  for  the  mass,  momentum, 
enthalpy,  and  Sauter-mean  diameter  variation,  Sp.  In  contrast  to  the  Lagrangian  formulation,  a 
single  averaged-size  particle/droplet  represents  the  dispersed  phase  which  makes  the  Eulerian 
analysis  computationally  more  efficient. 


&L. +ah.A +?2z 

dt  dt;  dri  dC, 


~  Fp,d  +  Fp,p  +  Fp,L 


(2.16) 


The  array  of  conserved  variables,  Qp,  and  the  flux  in  the  \  direction,  Ep,  are  given  by: 
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where  hp  is  the  particle  enthalpy  and  Sp  is  the  total  surface  area  of  the  particles  per  unit  volume.  The 
fluxes  Fp  and  Gp  are  defined  similarly.  The  source  terms  are  analogous  to  those  of  the  previous 
section. 


2.3  NUMERICS 

The  CRAFT  code  uses  a  finite-volume  upwind  formulation  with  an  implicit  (ADI) 
integration  scheme.  The  details  of  the  numerics  for  a  gas-phase  have  been  given  in  Ref.  4  while  the 
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fundamentals  for  extending  the  numerics  to  a  gas  liquid  system  have  been  summarized  in  Ref.  5. 
The  flux  at  cell  faces  are  using  a  Roe-TVD  flux  differencing  scheme  [6].  The  inclusion  of  the  gas- 
liquid  thermochemistry  described  above  into  the  upwind/implicit  methodology  is  complex  and 
requires  the  re-evaluation  of  the  flux  Jacobians,  as  well  as  the  eigenvalues  and  eigenvectors.  The 

Jacobian  of  the  flux  vectors  requires  the  evaluation  of  the  derivatives  such  as  etc.  To 

H  dPm 

evaluate  the  derivatives  of  pressure  we  require  the  functional  dependence  of  pressure  for  the  gas  and 
liquid  phases: 


Pg  fiPg’T’Qg)’  dPg  f[dpg,dT,d<Pg ) 

pL = /(Pi.r.O  -  Pj)).  dPL = /{dPodr.-d^)  (2'18) 

Inspection  of  Eqn.  (2.18)  reveals  that  while  the  pressure  differential  depends  on  (d<J>g),  4>g  itself  is 

not  a  variable  in  the  vector  Q.  Therefore  the  primary  obstacle  in  determining  the  pressure  derivative 
is  in  eliminating  the  dependence  of  the  pressure  differential  on  (d<j)g).  This  is  done  by  enforcing  the 
differential  form  of  Amagat's  law,  namely: 

dPg  =  dPL 

(2  19) 

d<j)g  =  G(dpL,dpg,dT ) 


Eqn.  (2.19),  is  substituted  back  into  Eqn.  (2.18)  to  obtain  the  pressure  differential  in  terms  of  the 
variables  solved  in  the  vector  Q  as  follows: 

dP  =  dPg  =  dPL  =  H(dpL,dpg,dT)  (2.20) 

The  details  of  the  flux  Jacobian  A  are  given  in  Ref.  7. 

The  acoustic  speed  of  the  two-phase  system  is  determined  from  the  eigenvalues  of  the  flux 
Jacobian  (see  Ref.  7  for  details)  and  is  given  by 


_j _ i  ti  i  iL 

p„cl  y  Eftc,2  Pi?L 
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(2.21) 


As  Eqn.  (2.21)  indicates,  the  acoustic  speed  in  a  two-phase  gas/liquid  mixture  behaves  differently 
than  in  either  a  pure  gas  or  a  pure  liquid.  For  two-phase  mixtures,  the  acoustic  speeds  of  the 
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individual  phases  do  not  combine  in  a  linear  fashion,  but  instead  exhibit  a  harmonic  relationship. 
The  variation  of  the  acoustic  speed  with  the  volumetric  fraction  of  the  gas  (for  a  given  pressure  and 
temperature)  has  a  "bathtub"  shape;  while  the  values  at  the  two  limits  are  the  respective  single  phase 
values,  the  acoustic  speed  drops  sharply  away  from  these  limits.  The  strong  dependence  of  the 
acoustic  speed  on  the  volumetric  fraction  puts  more  stringent  requirements  on  the  numerical 
algorithm  since  it  is  now  required  to  propagate  waves  with  the  correct  speed  in  mixtures  where  the 
two-phase  composition  is  varying  rapidly  (e.g.  across  a  gas/liquid  mixing  layer. 

2.4  EVAPORATION  MODELING 

Depending  on  the  droplet's  environment,  three  physical  models  have  been  identified  that 
describe  the  evaporation  process. 

2.4.1  Surface  Equilibrium,  Mass  Diffusion  Driven,  Corrected  for  Convection 

Under  most  circumstances,  the  duration  of  the  evaporation  process  is  much  greater  than  the 
gas-phase  time  scales,  and  the  vaporization  of  the  liquid  can  be  described  assuming  steady  state  in 
the  gas  phase  and  thermodynamic  phase-equilibrium  at  the  surface.  Unless  the  outside  pressure  is 
very  low  (very  high  release  altitude)  or  the  drops  are  very  small,  in  which  case  the  molecular  mean 
free  path  is  of  the  order  or  greater  than  the  drop  diameter,  the  rate  determining  process  is  the 
droplet  vapor  diffusion  in  the  gas  phase.  The  vaporization  rate  is  then  given  by  the  rate  of  mass 
diffusion  of  droplet  species  into  the  surrounding  gas.  This  rate  is  computed  by  assuming  that  the 
drop  is  immobile  with  corrections  applied  to  take  into  account  convection  due  to  the  relative 
velocity  between  the  drop  and  its  surrounding  [8].  Details  of  this  model  will  be  given  when  we 
described  the  SDROP  code  in  Section  5. 

2.4.2  Surface  Non-Equilibrium,  Kinetic  Rate 

For  very  low  surrounding  pressures  and  very  small  droplets,  the  assumption  of 
thermodynamic  equilibrium  at  the  surface  is  not  valid  and  the  vaporization  rate  is  computed  at  the 
drop  surface  as  the  difference  between  the  flux  of  molecules  of  drop  species  leaving  the  surface  to 
that  of  molecules  striking  the  surface.  This  kinetic  rate  is  given  by  the  Langmuir-Knudsen 
relationship  [9].  In  addition,  the  saturated  pressure  of  the  liquid  is  corrected  to  take  surface 
tension  into  account. 

2.4.3  Cavitation  Modeling  for  Flash  Vaporization 

In  cases  where  the  ambient  pressure  is  so  low  that  the  vapor  pressure  of  the  drop  species  at 
the  drop  temperature  is  greater  than  the  surrounding  pressure,  the  liquid  will  cavitate  and  bubbles 
will  form  inside  the  drop.  The  rate  of  vaporization  from  the  surface  is  computed  as  in  2.4.2  and  an 
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additional  equation  is  solved  for  the  birth  and  the  growth  of  gas  bubbles  inside  the  drop  [10].  Both 
the  mass  and  the  diameter  of  the  drop  are  adjusted  to  take  into  account  the  presence  of  the 
occlusions. 

2.5  COMPUTATIONAL  FRAMEWORK  AND  BREAKUP  MODELING 

The  overall  computational  framework  includes  three  separate,  coupled  integration  schemes. 
As  mentioned  in  the  Section  2.2.1,  the  motion  of  the  gas  and  the  liquid  is  computed  assuming  that 
the  gas/liquid  mixture  is  in  local  equilibrium  at  the  interface  for  which  we  apply  an  interface¬ 
capturing  Riemann  solver  to  the  NS  equations  for  the  mixture  interactions.  To  simulate  and  take 
into  account  the  effects  of  breakup,  the  gas/liquid  solver  is  fully  coupled  to  a  Lagrangian  and/or 
Eulerian  solver  for  the  droplet  phase.  In  this  section,  we  first  present  how  the  computational 
framework  is  organized  and  then  give  detailed  descriptions  of  the  solvers  and  the  physical  models 
they  employ. 

2.5.1  Description 

A  flow  chart  of  the  overall  computational  framework,  to  analyze  the  complete  problem 
schematized  in  Fig.  1.1,  is  shown  in  Fig.  2.1.  At  the  beginning  of  each  integration  time  step,  the 
code  computes  the  surface  tension  forces  and  checks  if  breakup  occurs  in  the  gas/liquid  interface 
cells.  In  cells  where  breakup  takes  place,  new  Lagrangian  parcels  are  generated.  The  coupling 
terms  between  the  discrete  (Lagrangian)  phase  and  the  gas/liquid  mixture  are  then  computed  and  the 
Lagrangian  parcels  are  updated.  If  the  number  of  Lagrangian  parcels  is  excessive,  Lagrangian 
parcels  can  be  transferred  to  the  Eulerian  particulate  continuum.  Then,  the  aerodynamic  drag  and 
heating  source  terms  resulting  from  the  interaction  of  the  gas  with  the  Eulerian  phase  are  updated, 
the  fluxes  for  the  Eulerian  particulates  are  computed,  and  the  Eulerian  particulate  variables  are 
updated.  Finally,  the  fluxes  for  the  gas/liquid  mixture  are  then  obtained  and  the  mixture  variables 
are  updated,  which  concludes  an  integration  cycle. 

2.5.2  Interface  Capturing  and  Surface  Effects 

As  explained  earlier,  the  fundamental  premise  of  our  breakup  work  is  that,  if  the  gas/liquid 
interface  can  be  sharply  captured,  then  surface  phenomena  can  be  simulated.  Therefore,  as  a  first 
step  towards  performing  numerical  simulations  of  liquid  breakup,  the  original  flux-limiting  (TVD) 
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methodology  of  Chakravarthy  and  Osher  [11]  was  redefined  in  terms  of  the  ratios  _  -  and 
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t  in  order  to  implement  varied  flux  limiters  and  test  their  effectiveness  in  limiting  the  numerical 
° .  1 

diffusion  at  the  interface.  Superbee  limiters  [12]  were  added  to  the  minmod  limiter  of  Ref.  1 1  and, 
for  the  simulations  of  Section  3.0,  the  Superbee  limiter  was  found  to  yield  the  sharpest  interfaces. 

Next,  the  CRAFT  code  was  updated  to  include  the  effects  of  surface  tension.  For  a  curved 
interface  between  two  fluids,  surface  tension  results  in  increased  pressure  on  the  convex  side  of  the 
surface  of  separation.  The  curvature-generated  pressure  gradient  can  couple  with  external 
excitations,  such  as  aerodynamic  forces,  to  amplify  surface  waves  and  cause  liquid  breakup. 
Surface  tension  may  be  an  important  parameter  in  controlling  the  deformation  and  breakup  of  liquid 
blobs,  and  should  be  represented  accurately  in  numerical  simulations  of  such  phenomena. 

The  relationship  between  surface  tension  and  the  stresses  at  the  interface  is  derived  by 
enforcing  mechanical  equilibrium  at  the  surface,  so  that  the  forces  on  both  sides  of  the  interface 
balance  each  other  [13].  At  the  interface  between  two  fluids: 


P2-a 


u 


(2.22) 


Eqn.  (2.26)  is  written  for  a  surface  which  is  convex  on  the  side  of  fluid  1,  with  P  representing 
pressure,  and  R2  the  principal  radii  of  curvature  of  the  surface,  a  the  surface  tension  coefficient, 
o'  the  viscous  stresses,  and  n-  the  components  of  the  surface  normal  in  the  three  directions.  For 
inviscid  liquids  with  constant  a,  Eqn.  (2.22)  reduces  to  Laplace’s  formula: 


Pi -P)=a\ 


V*1  R2J 


(2.23) 


In  numerical  methods  where  the  fluid  interface  is  tracked  explicitly,  Eqn.  (2.22)  is  simply  enforced 
as  a  boundary  condition  on  the  surface.  However,  when  the  interface  between  the  fluid  is  captured, 
not  fitted,  techniques  must  be  devised  to  include  the  effects  of  surface  phenomena  in  this  volume- 
based  description. 

The  CRAFT  code  has  been  updated  by  including  the  effects  of  surface  tension  using  the 
Continuum  Surface  Force  (CSF)  model  [14].  In  Brackbill  et  al.’s  formulation,  the  surface  tension 
force  is  recast  as  a  volume  force  which  is  added  to  the  momentum  equations.  The  magnitude  of  the 
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force  and  the  cells  to  which  it  applies  is  computed  based  on  the  gradient  of  a  function,  c,  defining 
the  location  of  the  different  fluids: 


F„(5)-a«:(x)^  (2.24) 

where  K  is  the  surface  curvature  and  c  the  jump  in  [c]  across  the  interface.  Ideally,  c  should  be  a 
smoothly  varying  function  such  as  a  distance  function,  however  it  can  also  be  taken  as  a  color 
function,  constant  in  each  fluid  and  changing  from  0  to  1  across  the  interface.  In  the  calculations 
presented  in  Section  3.0,  the  color  function  was  taken  as  the  mass  fraction  of  liquid,  or,  alternatively, 
the  volume  fraction  of  liquid. 

The  curvature,  K,  can  be  shown  to  be  equal  to: 
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where  the  vector,  n  =  Vc,  is  normal  to  the  liquid  surface,  pointing  towards  the  liquid.  Since  both 
the  magnitude  and  the  domain  of  application  of  Fsv  depend  on  K,  great  care  must  be  taken  when 
evaluating  n.  In  the  first  implementation,  V.n  was  computed  using  the  cell-centered  normals  and 
central  differencing.  However,  in  agreement  with  Brackbill  et  al.,  it  was  found  that  this 
implementation  resulted  in  an  increased  stencil  for  the  computation  of  Fsv,  in  effect  diffusing  the 
interface.  The  formulation  was  then  changed  to  compute  the  divergence  based  on  the  cell-vertex 
normals,  which  are  obtained  by  linear  interpolation.  It  was  also  found  that  the  partial  derivatives  of 
the  magnitude  of  the  normal,  n 
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in  order  to  be  consistent  with  V.«.  This  new  implementation  was  tested  on  cylindrical  and  spherical 
geometric  shapes  and  found  to  yield  accurate  curvatures.  After  the  interface  normals,  curvature,  and 
surface  tension  forces  are  computed,  the  forcing  term  in  Eqn.  (2.1)  is  expressed  as: 


S  -  [O’ Fsv,x ’Fsv,y> Fsv,z ,  0,  •  •  • ,  0,  •  •  • ,  0, 0] 


(2.28) 


The  implementation  of  the  CSF  model  was  tested  by  computing  the  deformation  of  an 
initially  cylindrical  pellet  of  liquid  propellant.  Results  for  the  liquid  mass  fraction,  mixture  velocity- 
vector,  and  pressure  field  are  shown  in  Fig.  2.2  at  three  different  times.  The  axis  of  revolution  is 
parallel  to  the  X  axis  and  cuts  are  presented  in  the  XY  plane.  Originally,  the  mixture  is  immobile 
and  both  the  temperature  and  pressure  are  uniform  throughout  the  computational  domain  so  that  the 
flow  and  the  pressure  gradient  that  develop  in  time  uniquely  result  from  the  effects  of  surface 
tension.  In  this  calculation,  the  initial  pressure  is  large  (100  atm)  so  that  the  liquid  is  slightly 
compressible.  In  order  to  exaggerate  and  observe  the  effects  of  the  surface  tension  force,  a  was 
artificially  augmented  accordingly.  The  surface  tension  coefficient  was  chosen  so  that  the 
equilibrium  pressure  in  the  liquid  would  be  about  1.1  times  that  in  the  gas. 

As  shown  by  the  velocity  vectors,  the  comers  are  set  in  motion  and  move  inward  toward  the 
X  axis,  while  the  top  of  the  drop  moves  in  the  outward  direction,  resulting  in  a  rounding-up  of  the 
pellet.  The  pellet  goes  through  several  oscillation  cycles,  during  which  it  successively  expands  and 
contracts  along  the  X  and  Y  axes,  slowly  converging  towards  a  spherical  shape  and  a  uniform 
pressure  inside  the  liquid.  Having  verified  the  implementation  of  the  CSF  model,  next  a  model  was 
developed  to  simulate  primary  breakup. 

2.5.3  Primary  breakup 

In  all  breakup  regimes,  the  primary  breakup  of  a  blob  of  liquid  is  preceded  by  the 
deformation  of  the  liquid  and  the  formation  of  surface  waves.  The  formation  of  drops  of  the  order 
of,  or  one  or  two  orders  of  magnitude  less,  than  the  bulk  diameter  can  be  resolved  by  computing  the 
merging  of  isosurfaces  of  liquid  mass  fraction,  and  this  may  be  enough  to  describe  breakup  in  the 
Rayleigh  regime,  where  drop  sizes  are  substantive.  However,  in  more  general  situations,  for  high¬ 
speed  flows,  the  byproducts  of  primary  breakup  can  be  several  order  of  magnitude  smaller  than  the 
bulk,  and  the  computational  grids  required  to  resolve  the  surface  waves  will  be  too  large  for  detailed 
breakup  calculations  to  be  conducted.  To  circumvent  this  difficulty,  models  are  used  to  compute  the 
rate  of  formation  of  droplets  and  to  characterize  their  size. 

Different  models  have  been  derived  to  simulate  primary  breakup  based  on  the  flow  and 
liquid  characteristics.  Reitz  [15]  has  obtained  correlations  for  the  rate  of  breakup  and  the  most 
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likely  drop  size,  dp.  Empirical  correlations  have  also  been  derived  for  the  breakup  time,  dp,  and  jet 
liquid  core  length  [16].  The  breakup-time  correlations  may  be  useful  for  well  defined  problems, 
such  as  the  interaction  between  a  liquid  drop  and  a  shock  wave.  However  they  require  one  to  be 
able  to  identify  the  time  zero  for  breakup,  which  limits  their  applicability.  The  correlations  of  Reitz 
[15]  have  proven  successful  in  simulations  of  liquid  injectors  [17],  however,  they  make  use  of  liquid 
length  scales  to  compute  the  We  and  Oh  numbers,  which  restricts  their  application  to  problems 
where  the  original  shape  of  the  liquid  is  well  known  and  a  length  scale  readily  identifiable.  This  is 
clearly  not  the  case  with  the  liquid  blob  breakup  problem  of  Fig.  1.1. 

For  more  general  problems  with  arbitrary  shapes  and  to  predict  the  onset  of  breakup,  a 
general  expression  for  the  rate  of  mass  breakup  at  the  surface  is  better  suited.  The  semi-empirical 
correlation  of  Mayer  was  implemented  in  the  CRAFT  code  as  formulated  in  the  Coaxial  Injection 
Combustion  Model  (CICM)  [18].  The  stripping  rate  takes  the  form: 


p1(p„f/2r)2  - 

mb  =  CA[  U  ^  ]3rcP,(A z) 

<*i/Pi 


and  the  drop  diameter  is  computed  as: 


i'.BjtdsiELf 

P  PgU2r 


(2.29) 


(2.30) 


The  empirical  constants,  CA  and  SA,  are  equal  to  0.08  and  120,  respectively,  as  advised  (p.  76)  in 
Ref.  18.  The  subscript  l  denotes  liquid  quantities,  g  denotes  gas  quantities,  p  is  the  viscosity  and  UT 
is  the  magnitude  of  the  relative  velocity  between  the  two  phases.  The  CICM  was  originally 
developed  to  simulate  the  breakup  of  cylindrical  jets,  for  which  D,  represents  the  jet  diameter  and 
JtD,(Az),  the  surface  area  of  the  liquid  in  a  computational  cell  of  length  A z.  For  our  applications, 
Eqn.  (2.29)  was  generalized  by  reformulating  the  stripping  rate  as  the  mass  rate  of  breakup  per 
surface  area  of  the  liquid,  5,: 


mh  MpsU2r)  i 


(2.31) 


Knowing  how  to  compute  and  d  ,  the  outstanding  difficulty  in  applying  the  surface 


breakup  model  is  to  decide  where  and  when  breakup  occurs.  To  that  end,  and  to  simulate  the 
physical  mechanism  responsible  for  breakup  which  is  that  surface  waves  grow  until  they  become 
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unstable  and  break  off,  surface  breakup  was  modeled  as  a  rate  process,  based  on  the  local, 
computed,  properties  of  the  captured  interface.  The  strength  of  this  approach  is  that  the  breakup  is 
modeled  based  on  local,  dynamic  properties,  not  on  a-priori  or  global  considerations.  Therefore, 
provided  a  fundamental  correlation  is  used  to  compute  the  breakup  rate,  the  physics  of  breakup  will 
be  represented  accurately. 

The  present  model  is  obtained  by  analogy  with  chemical  kinetic  modeling,  where  a  reaction 
mechanism  is  divided  into  an  induction  and  a  reaction  period  whose  evolution  is  followed  by 
computing  the  change  in  the  fraction  of  elapsed  induction  and  reaction  times.  In  the  present  case, 
we  track  the  mass  available  for  breakup,  and  breakup  occurs  when  enough  mass  has  been 

accumulated.  Defining  the  fraction  of  broken  mass,  /  =  — ,  where  mb  is  the  mass  available  for 

md 

breakup  and  md  the  mass  of  a  drop,  we  compute: 


dt  md 


(2.32) 


The  derivative  in  Eqn.  (2.32)  is  a  particle  derivative  and/is  convected  with  the  liquid  so  that: 


+  m.V/  =  —  (2.33) 

at  md 


Taking  into  account  the  breakup,  the  conservation  of  mass  equation  becomes: 

+  V.(wp)  =  m 
dt 


(2.34) 


where  the  rate  of  mass  change  of  the  mixture,  m ,  in  Eqn.  (2.34)  is  computed  as  the  total  mass 
breaking  up  in  a  computational  cell  during  an  integration  time  step.  Combining  Eqns.  (2.33,2.34),  a 
conservation  equation  for/ is  obtained: 

^  af) +  =  P~  +  fi1  (2.35) 

at  md 

Eqn.  (2.35)  was  implemented  into  the  CRAFT  code,  as  described  in  Section  2.2,  with  the  resulting 
source  term: 
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(2.36) 


S  =  A tVc  m,mud,mvd,mwd,mhd,0,"-,m,p—  +  fih 
L  md 

Eqn.  (2.36)  implies  that,  during  a  timestep  At,  where  breakup  occurs,  a  mass  mx  At  of  liquid  is 

removed  from  the  mixture,  along  with  its  momentum  and  enthalpy,  and  p—  +  fin  less  liquid  mass 

md 

is  available  for  breakup.  With  this  new  model,  at  each  step,  the  fraction  of  broken  mass  is 
computed  in  each  cell  and: 

2  if /<  1,  no  breakup, /increases  according  to  Eqn.  (2.32) 

3  if />  1,  form  Nd  =  int (/)  drops  of  nominal  diameter  given  by  Eqn.  (2.30). 

The  final  step  in  implementing  our  model  consists  of  computing  the  dynamic,  local  parameters, 
namely,  pg,  p,,  Ur,  and  5,,  needed  in  Eqns.  (2.30)  and  (2.31).  In  keeping  with  our  philosophy  of 
capturing  the  gas/liquid  discontinuity  and  distributing  surface  phenomena  over  the  numerical  width 
of  the  interface,  stripping  is  applied  to  all  the  computational  cells  in  which  e  <  c  <  1-e,  where  £  is  a 
small  adjustable  number  (typ.  1%). 

Particular  care  must  be  taken  when  computing  5,  in  order  to  avoid  generating  exceedingly 
large  breakup  fluxes.  Consider,  for  example,  the  simulation  of  a  cylindrical  jet  of  radius  R. 
Analytically,  the  surface  area  of  liquid  is  tcR(Az).  In  the  calculations,  however,  the  interface  will 
spread  from  (R-Ar)  to  (R+Ar)  due  to  numerical  diffusion.  If  the  discontinuity  is  spread  over,  say,  3 
cells,  the  total,  computed  surface  area  will  then  be  n(R-Ar+R+R+Ar)(Az),  or  three  times  the  exact 
surface  area.  Therefore,  the  surface  area  in  each  cell  must  be  weighed  to  ensure  that,  in  the  limit  of 
an  infinitely  small  thickness,  5,  converges  to  the  exact  value.  This  is  accomplished  in  the  present 
methodology  by  analogy  with  the  CSF  model. 

By  definition,  the  surface  tension  force  per  unit  area  of  the  liquid  is  equal  to: 

Fsa(x)  =  aK(x)n(x)  (2.37) 

where  the  magnitude  of  n(x)  is  equal  to  the  surface  area.  As  described  in  Section  2.5.2,  in  the  CSF 
model,  is  recast  as  the  volume  force  Fsv.  Therefore,  to  insure  the  proper  weighing  we  enforce: 

FsaS[  =  FsvVc  (2.38) 

where  Vc  is  the  volume  of  a  cell  and  the  liquid  surface  area  in  each  cell  is  then  obtained  as: 
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(2.39) 


Implementation  of  the  dynamic  breakup  framework  was  completed  by  developing  a 
methodology  to  compute,  pg,  pv  and  UT.  Inherent  to  the  interface  capturing  methodology  is  the  fact 
that  the  gas/liquid  discontinuity  is  spread  over  a  few  cells  containing  both  gas  and  liquid  in 
equilibrium.  Therefore,  in  each  interface  cell,  a  nearest,  appropriate,  gas  cell  must  be  found  in  order 
to  compute  pg,  and  UT. 

In  the  present  framework,  we  assume  that  gas  cells  whose  properties  dictate  the  rate  of 
breakup  are  situated  nearest  and  directly  above  the  liquid  surface,  in  the  direction  of  n.  For  each 
interface  cell,  the  nearest  gas  cell  is  found  by  forming  the  cross-product  of  n  and  the  cell-centers 

vector  for  each  neighbor  cell,  CcQ,  as  shown  in  Fig.  2.3.  The  search  can  be  extended  to  several 
cells  around  the  breakup  cell,  and  the  chosen  gas  cell  is  that  for  which  the  angle  between  n  and 
CcCij  is  the  smallest  and  so  that  c  <  e.  The  latter  requirement  is  necessary  to  prevent  breakup  in 

areas  where  mixed  cells  may  be  surrounded  by  liquid  cells,  due  to  extreme  convolution  of  the 
interface,  as  happens  in  the  results  to  be  shown  in  Section  3.2.  The  same  operations  are  repeated  to 
find  the  pure  liquid  cell.  Once  the  gas  and  liquid  cells  are  found,  and  pg,  and  p,  are  computed,  Ur  is 

obtained  by  forming  the  cross  product  of  the  relative  velocity  vector  Ug  -  U[  and  n.  Alternatively, 

Ut  has  also  been  obtained  by  forming  the  cross  product  of  Ug  -  Ut  and  the  shear  stress  vector. 

In  every  interface  cell  where  breakup  occurs,  the  number  of  drops,  Nd,  formed  during  one 
integration  step  is  obtained  by  dividing  the  total  breakup  mass  by  the  mass  of  one  drop  of  nominal 
diameter.  A  new  computational  parcel  representing  Nd  physical  drops  is  then  initializfvl  at  the  cell 
center,  with  the  mixture  velocity  and  temperature  at  the  time  of  breakup.  In  future  work,  different 
initial  conditions  will  be  explored,  for  example  including  a  random  ejection  velocity  as  in  Amsden 
[19]. 
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Fig.  2.1.  Schematic  of  an  integration  time  step  for  complete  problem  shown  in  Fig.  1.1. 
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Fig.  2.2.  Deformation  of  a  cylinder  of  liquid  propellant.  Contours  of  liquid  isopleths  (top), 
velocity  vectors  (middle),  and  contours  of  pressure  (bottom). 


Fig.  2.3.  Schematic  of  the  computational  grid  in  two  dimensions 
showing  the  gas/liquid  interface. 
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3.0  CRAFT  CODE  BULK  LIQUID  FLYOUT  BREAKUP  STUDY 


3.1  OVERVIEW 

The  study  to  be  described  exhibits  the  ability  of  the  CRAFT  Navier-Stokes  code,  with 
enhanced  gas/liquid  and  new  breakup  capabilities,  to  analyze  the  deformation  and  breakup  of  a 
liquid  blob,  travelling  at  high  speed.  The  blob  has  an  initial  cylindrical  shape  with  the  flat  faces  at 
the  front  and  rear.  Conditions  and  size  correspond  to  earlier  experiments  performed  at  AEDC.  The 
blob  is  released  with  an  initial  velocity  of  700  m/s  (Mach  number  =  2)  and  has  an  initial  diameter 
and  length  of  50  cm.  The  grid  utilized  moves  with  the  flying  blob  in  accordance  with  the  velocity  of 
the  grid  point  at  its  center  of  gravity  (initially).  The  working  fluid  for  this  case  was  a  conventional 
liquid  propellant. 

The  complexities  of  this  flow  problem  can  be  gleaned  from  the  instantaneous  solution 
snapshot  at  .9  ms  (Fig.  3.1).  The  shape  has  deformed  to  have  flattened  teardrop  characteristics,  a 
blunt  body  shock  sits  in  front  of  the  decelerating  blob,  a  zone  of  separated  flow  is  evident 
downstream  of  the  shoulder,  and  the  disturbed  aerodynamic  flow  enters  the  wake  behind  the  blob 
and  is  heated  by  the  nose  shock  entropy  rise.  The  flow  within  the  blob  has  extremely  complex 
characteristics  due  to  the  pressure  variations  produced  by  the  aerodynamic  interactions. 

Two  sets  of  calculations  were  performed,  the  first  without  aerodynamic  breakup  and  the 
second  with  aerodynamic  breakup.  With  aerodynamic  breakup,  mass  is  released  in  accordance  with 
an  extension  of  spray-based  correlations  for  droplet  sizes  smaller  than  those  resolvable  by  the 
computational  grid  utilized.  It  would  be  feasible,  but  not  practical,  to  directly  simulate  the  formation 
of  droplets  at  the  interface,  which  would  require  grid  scales  substantially  smaller  than  the  scales  of 
the  droplets  formed.  As  will  be  seen,  droplet  globs  which  can  be  resolved  by  the  grid  do  form  at  the 
surface  and  break-off.  Smaller  droplets  leave  the  gas/liquid  interface  in  accordance  with  mass  rate 
and  droplet  size  relations  that  have  dependencies  on  liquid  properties  (viscosity  and  surface  tension) 
and  on  dynamic  interface  conditions  (gas/liquid  relative  velocity  and  density  ratios).  The 
implementation  of  these  relations  for  a  dynamic/deforming  blob  is  extremely  complex  and  we  refer 
the  reader  to  Section  2  for  details. 
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Fig.  3.1.  Cylindrical  blob  flyout.  Details  of  the  flowfield  at  t  =  0.9  [ms].  Velocity  vectors 
superimposed  on  contours  of  liquid  mass  fraction  (black)  and  pressure  contours  (color). 


In  the  work  described,  we  divide  the  problem  in  terms  of  resolvable  and  non-resolvable 
scales.  The  motion  and  thermodynamic  states  of  the  bulk  liquid  are  computed  directly  based  on 
the  Navier-Stokes  (NS)  equations,  extended  to  gas/liquid  mixtures  in  local  equilibrium  at  the 
interface.  Liquid  primary  breakup  is  modeled  by  forming  droplets  at  the  interface  between  the 
liquid  and  the  surrounding  gas  based  on  surface  characteristics,  and  empirical  and  analytical 
models.  Once  formed  the  motion,  energy  and  conservation  of  mass  in  either  a  Lagrangian  or 
Eulerian  manner.  The  deformation  of  the  drops  and  secondary  breakup  can  be  taken  into  account 
using  analytical  and  empirical  models,  as  will  be  described  in  Section  5. 


-25- 


3.2  BULK  LIQUID  FLYOUT  WITHOUT  BREAKUP 

Initially,  a  cylindrical  pellet  is  impulsively  started  at  M  =  2  from  right  to  left,  in  still  air  (Fig. 
3.2).  The  grid  consists  of  240  cells  in  the  X  direction  and  180  cells  in  the  Y  direction,  with  a 
200x140  cells  fine-grid  region  ( Ax  =  Ay  =5  mm)  around  the  blob  surrounded  by  a  coarser-grid 
region  where  the  cell  size  increases  progressively  with  distance  from  the  fine  grid  according  to  a 
geometric  progression.  In  the  following  figures,  the  axis  of  symmetry  is  horizontal  and  centered  on 
Y  =0.0.  The  grid  moves  with  the  velocity  of  the  pellets  initial  center  of  gravity. 


Fig.  3.2.  Free  flight  of  a  cylindrical  pellet  of  fluid.  Contours  of  gas  mass  fraction. 


-26- 


As  a  result  of  gas/liquid  interaction,  pressure  builds  up  in  front  of  the  drop  and  the  pellet 
starts  to  deform,  as  shown  at  t  =  0.4  ms.  The  front  shoulder  of  the  blob  starts  rolling  back  and 
liquid  is  pulled  out  of  the  back.  Later  in  time,  at  t  =  2.4  pm,  the  liquid  carried  from  the  shoulder  and 
the  back  reflect  off  of  the  centerline  and  spread  outward,  away  from  the  blob,  and  inwards,  towards 
the  back  side  of  the  pellet.  The  origin  of  these  complex  flow  features  can  be  elucidated  by 
simultaneously  looking  at  contours  of  pressure,  gas  mass  fraction,  and  at  the  velocity  vector,  as  in 
Figure  3.2. 

A  bow  shock  reflects  off  the  front  (left)  side  of  the  blob.  Air  is  entrained  behind  the  bow 
shock  and  in  the  wake  of  the  blob,  where  expansion  waves  accelerate  the  flow  from  right  to  left.  As 
the  bow  shock  is  formed,  a  planar  shock  is  transmitted  in  the  liquid,  increasing  the  pressure  in  the 
blob  and  causing  the  liquid  to  expand  laterally  (vertically),  as  seen  at  t  =  0.9  ms.  While  the  front 
half  of  the  blob  expands  laterally,  the  back  half  spreads  down  due  the  high  shear  and  the  expansion 
waves  in  the  back  comer  of  the  drop.  This  combination  of  expansion  and  shrinkage  give  the  blob 
the  shape  displayed  at  t  =  2.0  [ms],  with  a  flattened  disk-like  front  surface  connected  to  an 
elongated  horse-shoe  shaped  back  surface  through  an  inclined  interface.  The  lateral  expansion  of 
the  liquid  at  the  front  of  the  blob  causes  the  gas  to  move  vertically  and  sideways.  Thinking  in 
relative  coordinates,  with  the  reference  frame  attached  to  the  blob,  the  laterally  expanding  liquid  acts 
like  an  airfoil,  accelerating  the  gas  around  its  top  surface  and  causing  it  to  expand.  A  recirculation 
region  is  created  at  the  top  of  the  blob  and  gas  is  entrained  in  the  liquid.  Mushroom  structures 
characteristic  of  Rayleigh-Taylor  instabilities  are  visible  at  the  top  and  the  back  of  the  blob,  and 
distinct  regions  containing  only  gas  become  entrapped  in  the  liquid. 

This  analysis  shows  that  the  bulk  deformation  takes  place  through  a  combination  of 
complex  interactions  of  compression  waves  at  the  front  of  the  blob,  and  expansion  waves  and  shear 
at  the  back  of  the  blob.  As  a  result  of  these  interactions,  both  at  the  front  and  back  of  the  blob, 
lower-density  fluid  (gas)  is  accelerated  into  higher-density  fluid  (liquid),  raising  the  potential  for 
Rayleigh-Taylor  instabilities  to  set  in  and  act  as  a  driver  for  blob  breakup.  It  is  interesting  to  note 
that  Joseph  et  al.  [20]  have  already  alluded  to  such  a  possibility  based  on  their  experimental  results. 
These  results  indicate  that,  in  absence  of  surface  primary  breakup  modeling,  the  large  mass  of  fluid 
which  rolls  up  above  the  blob  would  tend  to  break  away  from  the  blob.  Next,  the  same  problem 
was  repeated  including  the  effects  of  primary  breakup. 

3.3  BULK  LIQUID  FLYOUT  WITH  BREAKUP 

As  can  be  seen  from  the  contours  of  gas  mass  fraction  in  Figure  3.3,  the  changes  with 
breakup  are  very  drastic.  The  large  roll-ups  observed  in  Figure  3.2  are  now  absent  and  vigorous 
mass  shedding  takes  place  at  the  shoulder  and  the  back  of  the  pellet,  where  the  shear  is  highest. 
The  computations  were  repeated,  this  time  including  two-way  coupling  with  the  broken-off  drops. 
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In  these  results,  the  rate  of  mass  breakup  is  very  high  and  a  large  number  of  computational  parcels 
are  generated  at  each  step.  Therefore,  all  the  Lagrangian  parcels  are  transferred  into  the  Eulerian 
continuum  at  each  time  step.  Figure  3.4  shows  the  Eulerian  cloud  of  particulates  which  forms  a 
trail  behind  the  pellet.  Very  high  particle  densities  are  observed  near  the  drop  shoulder  and  the  back 
of  the  drop,  in  qualitative  accord  with  photographs  of  shock  wave  /  drop  interactions  in  shock  tubes 
[20].  The  drop  size  distributions  are  also  in  qualitative  accord  with  experimental  observations  of 
drop  shattering.  The  smaller  drops,  which  equilibrate  fastest  with  the  slower  moving  gas,  are  found 
on  the  outer  edges  of  the  cloud  and  the  larger  drops  are  close  to  the  pellet.  In  the  future,  boundary 
conditions  inside  the  computational  domain  will  be  implemented  to  prevent  the  Eulerian  drops  from 
diffusing  inside  the  liquid,  and  the  effects  of  secondary  breakup  will  be  evaluated. 

The  shapes  predicted  appear  to  be  in  qualitative  accord  with  those  exhibited  in  laboratory 
studies  at  the  University  of  Minnesota  [20],  and  the  relatively  fast  rate  of  breakup  under  these 
conditions  appears  to  be  in  nominal  accord  with  gross  observations  from  the  AEDC  tests.  In 
particular,  Fig.  3.5  shows  the  rate  of  mass  breakup  with  time  for  this  case,  with  all  mass  broken  up 
in  6  ms.  Figure  3.6  shows  the  blob's  position  vs.  time  for  this  case  (and  for  the  blob  with  no 
breakup)  which  would  indicate  that  the  flyout  distance  for  complete  breakup  is  about  3  meters. 
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Fig.  3.4  :  Free  flight  of  a  cylindrical  pellet  of  fluid.  Breakup  with  Eulerian  tracking  of  the  particulate 
continuum. Contours  of  gas  mass  fraction  (black)  superimposed  on  top  of  contours  of  particulate 
cloud  density  (color,  left)  and  contours  of  cloud  SMD  (right,  color). 
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4.0  DROPLET  CLOUD  NEUTRALIZATION  STUDIES 


Exploratory  studies  were  performed  to  examine  the  response  of  a  cloud  of  droplets  to 
interactions  with  a  blast  that  is  set  off  in  their  path  of  propagation.  Droplet  clouds  having  uniform, 
spherical  droplet  sizes  of  100  pm  and  10pm  moving  with  an  initial  velocity  of  700  m/sec  were 
considered.  A  stationary  point  blast  was  set  off  just  prior  to  the  clouds  arrival.  Droplet 
decomposition  was  represented  by  a  "plausible"  kinetic  mechanism  requiring  a  surface  temperature 
of  400°K  to  become  activated.  Figures  4. la-4. Id  exhibit  contours  of  pressures  and  particle 
temperatures,  for  times  initiating  just  prior  to  the  cloud/blast  interaction  and  terminating  when  the 
droplets  have  propagated  mostly  through  the  blast.  As  can  be  gleaned  from  these  contours,  the 
blast  is  quite  effective  in  neutralizing  the  10pm  droplets  but  is  not  effective  in  neutralizing  the 
100pm  droplets.  Since  heat  transfer  is  proportional  to  surface  area,  the  smaller  droplets  can  absorb 
the  heat  produced  by  the  blast  at  a  much  faster  rate.  Timing  is  also  critical  since  peak  temperatures 
generated  by  a  simple  point  blast  will  rapidly  decay  (Figure  4.2)  and  effective  neutralization 
requires  exposure  to  high  temperature  for  a  time  which  suffices  for  surface  temperatures  to  reach  a 
critical  value,  for  the  kinetic  mechanism  implemented  in  this  study.  Figure  4.3  exhibits  the  rate  of 
neutralization  of  the  10pm  droplets  in  this  scenario  and  contrasts  it  with  the  rate  for  the  100  pm 
droplet,  which  does  not  decompose. 

There  are  a  number  of  important  points  to  be  gleaned  from  these  exploratory  studies, 

namely: 

(1)  first-principles  models  can  support  the  assessment  of  neutralization  concepts; 

(2)  droplet  size  needs  to  be  well  characterized  with  the  primary  breakup  work  described  in 
Section  2  providing  a  path  towards  simulating  post-hit  size  distributions; 

(3)  droplet  decomposition  mechanisms  (for  vaporization/combustion)  need  to  be  kinetic  based 
and  accurately  modeled  since  the  surrounding  environment  is  highly  dynamic;  and 

(4)  details  of  how  the  neutralization  is  performed  (timing,  charge  location,  etc.)  are  also  very 
important  since  one  needs  to  "capture"  the  droplets  in  a  sustained  hot  environment  and  be 
assured  that  the  blast  does  not  displace  them  from  the  hottest  zones. 

Our  current  Air  Force  work  in  simulating  bunker  neutralization  scenarios  [21]  has 
analogous  challenges.  In  such  scenarios:  canisters  filled  with  C/B  agents  are  punctured;  escaping 
agent  (e.g.  a  liquid  jet)  breaks-up  into  a  droplet  cloud;  sustained  thermal  energy  from  chemical 
mechanisms  is  required  for  neutralization;  and,  care  is  required  that  blast  pressurization  does  not 
vent  the  agent  outside  the  bunker  prior  to  neutralization.  Varied  studies  are  described  in  Ref.  [21] 
that  simulate  bunker  neutralization  scenarios  using  the  CRAFT  code,  for  which  data  has  been 
obtained  at  Southwest  Research  Institute  and  Eglin  Air  Force  Base. 
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Fig.  4.1a.  Blast  /  particulate  cloud  interaction,  t  =  3.5  [ms]. 
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Fig.  4.1b.  Blast  /  particulate  cloud  interaction,  t  =  4.0  [ms] 
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Fig.  4.1c.  Blast  /  particulate  cloud  interaction,  t  =  4.3  [ms]. 
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Fig.  4.1d.  Blast  /  particulate  cloud  interaction,  t  =  4.6  [ms]. 
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5.0  SDROP  METHODOLOGY  AND  PARTICLE  FATE  IMPLICATIONS 


5.1  OVERVIEW 

In  high  altitude  post-hit  scenarios,  the  fate  of  individual  droplets  is  predicted  in  Systems 
Models  such  as  PEGEM  and  VLSTRACK,  by  simple  trajectory  models  which  use  drag  and  heat 
transfer  correlations,  supplemented  by  rudimentary  vaporization  and  aero-breakup  models. 
Concerns  were  expressed  regarding  the  accuracy  of  such  models,  with  predictions  made  by 
different  groups  leading  to  very  different  estimates  of  particle  fate.  Our  role  in  this  endeavor  was  to 
develop  a  higher-fidelity  engineering  model  that  could  provide  improved  predictions  and  serve  to 
indicate  what  parameters  in  Systems  Models  needed  to  be  upgraded.  The  resultant  model,  SDROP, 
implemented  the  Lagrangian  particle  methodology  in  the  CRAFT  code,  coupled  with  3  DOF 
trajectory  equations  and  atmospheric  tables. 

A  Review  of  Systems  Model  methodology  indicated  potential  deficiencies  regarding: 
drag/heat  transfer  relations,  the  assumption  of  uniform  droplet  temperature,  the  vaporization  models 
implemented,  and  the  aerobreakup  methodology  utilized.  In  SDROP,  advanced  models  for  the 
above  were  utilized  as  summarized  in  Table  5.1. 


Table  5.1  -  Droplet  Fate  Model  Improvements  in  SDROP 


Parameter  Modeled 

Approach  Utilized 

Drag/Heat  Transfer 

Rocket  exhaust  particulates  containing  lag 

Mach/Reynolds  number  corrections  (Ref.  3) 

Droplet  Temperature 

ODE  forT(r,t)  and  analytical  approximation 

Surface  Vaporization 

Equilibrium  (Spalding,  Ref.  22)  and  nonequilibrium 
models  with  coupled  heat  transfer 

Aerodynamic  Breakup 

Droplet  deformation  model  of  Liu  and  Reitz  (Ref.  23) 

Varied  studies  that  we  have  performed,  to  be  described  below,  have  indicated  the  importance 
of  using  more  advanced  approaches  for  predicting  droplet  fate.  We  have  identified  sensitivities  to 
both  aerodynamic  breakup  modeling  and  to  use  of  variable  temperature.  The  simple  Systems 
Model  approach  does  not  consider  shape  deformation  and  uses  a  simplistic  Weber  number*  (We) 


2 

dynamic  force  2rpPg  vp  -  v„  ,  .  .  .  ,  .  i 

=  — — - = - where  r  is  droplet  radius,  p  is  gas  density,  |  v  -v 


*  We  = 

surface  tension  force 
is  relative  velocity  and  as  is  surface  tension  as  the  sole  basis  for  breakup. 


p  g 1 
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criterion.  In  this  criterion,  when  the  critical  stability  value  (We~12)  is  reached,  aerodynamic 
breakup  occurs.  In  the  advanced  approach  utilized,  the  droplet  first  deforms  in  accordance  with  the 
deformation  rate  equations  of  Ref.  23  which  assume  an  ellipsoidal  shape.  Breakup  occurs  when 
deformation  produces  a  major/minor  axis  deformation  ratio  that  exceeds  2.  As  we  will  discuss 
below,  for  a  1  cm  VX  droplet  released  at  90  km,  the  critical  Weber  number  is  reached  within  10 
seconds  at  an  altitude  of  68  km  and  the  prediction  indicates  no  droplets  reaching  earth.  In  contrast, 
the  deforming  droplets  maintain  stability  due  to  changes  in  the  aerodynamic  drag  produced  by 
deformation  and  all  droplets  reach  earth.  Variable  temperature  within  the  droplet  also  produces 
changes  which  have  significant  impact.  Droplets  do  not  heat/cool  uniformly  and  the  surface 
temperature  effects  transport  coefficients.  In  a  vaporizing  case,  the  effects  would  be  dramatic.  For 
a  non-volatile  simulant  such  as  VX,  where  vaporization  plays  a  negligible  role  in  droplet  fate,  the 
effect  is  relatively  small. 

5.2  SDROP  DESCRIPTION 

SDROP  is  a  stand-alone  Lagrangian  computer  code  formulated  to  calculate  the  fate  of 
liquid  droplets  released  in  the  atmosphere  for  given  initial  conditions.  The  particular  impetus  of 
this  work,  made  possible  by  the  modular  nature  of  the  code,  is  to  provide  a  tool  to  evaluate  varied 
physical  models  and  the  relative  importance  of  the  physical  phenomena  involved,  such  as 
evaporation  and  breakup.  Given  the  initial  release  altitude,  diameter,  temperature,  and  velocity 
vector  of  a  drop,  SDROP  computes  its  three-dimensional  trajectory,  taking  into  account  the 
changes  in  atmospheric  properties,  heat  and  mass  transfer  between  the  drop  and  its  surroundings, 
and  drop  deformation  and  breakup  due  to  aerodynamic  forces.  Details  of  varied  components  of 
SDROP  are  provided  in  the  subsections  that  follow. 

5.2.1  Drop  Trajectory 

The  change  in  position  and  velocity  of  the  drop  is  computed  in  earth-fixed  coordinates. 


d2x  1  , 

^4  = — Fx+2Vvco  +  xco2 
dt2  mp  y 


(5.1) 


-  2  Vx(o  +  yo)2 


(5.2) 


d~Z=±_F 

dt2  mp  z 


(5.3) 
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co  is  the  rate  of  rotation  of  the  earth.  In  the  present  implementation  of  SDROP,  the  force  on  the 
droplet  is  computed  for  the  following  set  of  assumptions: 

1 .  Spherical  droplets  behave  like  solid  spheres 

2.  Droplets  do  not  rotate 

3.  Atmospheric  turbulence  neglected 

4.  Axisymmetric  uniform  flow  around  droplets 

4a.  Buoyancy  force  and  aerodynamic  drag  only,  no  lift 

4b.  Neglect  Basset  history  term  and  added  p  mass  terms,  which  is  standard  procedure 
when  — »1 

Ps 

and  F  is  then  obtained  from: 


F  = 


\psAfCD  p-Vp-V) 


+  mp 


(5.4) 


5.2.2  Representation  of  the  Atmosphere 

In  the  results  presented  in  Section  5.3,  data  on  wind  velocity  was  assumed  to  be  zero. 
(U=0).  However,  approximate  curve  fits  taken  from  the  U.S.  Standard  Atmosphere,  1962  [24]  are 
used  to  compute  the  variation  of  atmospheric  temperature,  density,  pressure,  molecular  weight,  and 
viscosity  with  geometric  altitude.  The  approximate  formulae  are  detailed  in  Appendix  A. 
Implementing  extended  or  additional  curve  fits  in  SDROP  would  be  straightforward  since  they  are 
written  in  the  form  of  Fortran  functions.  It  is  important  to  remember  that  “...the  U.S.  Standard 
Atmosphere,  1962  is  an  idealized  middle-latitude  (approximately  45)  year-round  mean  over  the 
range  of  solar  activity  between  minima  and  maxima.”  Therefore,  our  current  representation  of  the 
atmosphere  does  not  take  into  account  temporal  and  spatial  variations,  which  can  be  significant,  and 
upgraded  atmospheric  data  would  be  required  in  order  to  conduct  accurate  drop  fate  calculations  for 
realistic  release  scenarios.  However,  as  stated  in  the  overview,  SDROP  is  intended  to  conduct 
parametric  and  modeling  studies,  and  the  approximate  formulae  of  Appendix  A  provide  a 
sufficiently  realistic  description  of  the  atmosphere  to  achieve  that  goal. 

The  expression  for  the  acceleration  of  gravity  given  by  Ref.  24  includes  the  centrifugal 
acceleration  and  it  was  modified  since,  in  SDROP,  the  equation  of  motions  are  already  solved  in 
earth-fixed  coordinates 
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\s\  =  |VOc|  =  ~~\l  +  /(^)2(3sinyr2  -  i)||l  +  ^(^)8(27siny/-cosi//')J  (5.5) 


5.2.3  Aerodynamic  Drag 

The  drag  coefficient  was  evaluated  with  Henderson’s  formula  (see  Ref.  3),  where  Mp  and 
Rep  are  the  lag  Mach  number  and  Reynolds  numbers  (based  on  Up  -  Ug)  which  correct  the 
limiting  value  of  CD  of  24/Re  for  Mp  and  Rep  going  to  zero  (Stokes  flow). 


Mp  <  1.0 


24 


Mp  >  1.75  :  CD  = 


1.0  <  Mp  <1.75 


Q  — _ 

D  Rep+  s(4.33  +  hx  exp[-.247  Rep/  s]) 

+{h1  +.\M2p  +  .2M|)exp[-Mp/2^/Re^] 

+.6^1  -  exp[-Mp  /Rep])s 

.9  +  .34/  M2p  + 1 .86(Mp  /  Rep)^r2  +  2  /  S2  + 1 .058 /S(TP  -MS4 


\  +  l.S6(Mp/Repf2 

c  =  c  I  K-‘)  (r  _r  1 

D  Ldmp= 1  1.75-1  [Cdmp=i.15  Cdmp=i) 


(5.6) 


.  ,  3.65- 1.53(2}, /I) 

where :  h  = - } 

1+.353 (TP/T) 

4.5 +.38(.03Ren+  ,48ReP1/2) 

h2= - i - - - yp: — - 

1  l  +  .OSRep+^Rep172 


Ml 

2)  Mp 


5  =  W  M, 


While  several  correlations  are  implemented  in  SDROP,  the  Cd  and  Nu  correlations  that  are 
currently  used  were  derived  for  rocket-plume  applications  [3].  The  advantage  of  these  correlations 
compared  to  the  ones  often  used  for  conventional  applications  is  that  they  were  derived  over  a 
wider  range  of  lag  Reynolds  and  Mach  numbers  and  include  corrections  for  compressibility 
effects.  The  Cd  correlations  have  also  been  shown  to  be  valid  up  to  the  free-molecular  flow 
regime.  When  computing  the  drag  and  heat  transfer  to  the  drop,  corrections  are  also  included  to 
take  into  account  the  blowing  due  to  vaporization.  The  heat  transfer  coefficient  is  further  modified 
using  film  boiling  correlations  for  flash  vaporization.  Finally,  the  heat  transfer  is  based  on  the 
recovery  temperature  in  the  boundary  layer  around  the  droplet  and,  when  flying  supersonically,  the 
gas  properties  around  the  drop  are  computed  behind  the  shock  wave  generated  by  the  drop.  The 
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total  drag  on  the  droplet  is  computed  as  a  function  of  Cd  and  the  frontal  area,  both  of  which 
depend  on  the  shape  of  the  drop.  In  absence  of  deformation,  sphere  values  are  used.  When 
deformation  is  computed,  the  frontal  area  and  Cd  are  computed  using,  respectively,  the  true  frontal 
area  of  the  drop  and  a  linear  interpolation  between  a  disk  and  a  sphere. 

5.2.4  Mass  and  Heat  Transfer 

The  equations  governing  the  mass  and  heat  transfer  between  the  drops  and  the  surrounding 
atmosphere  are  presendy  valid  for  the  following  assumptions: 

1 .  Monocomponent  droplet  -  solve  only  for  the  droplet  mass 

2.  No  chemical  reactions  in  the  liquid  and  the  gas,  mass  transfer  occurs  only  as  result 
of  vaporization,  condensation  or  breakup 

3.  Liquid  motion  inside  the  drop  is  neglected 

4.  Two  models  are  available  for  the  particle  temperature: 

4a.  Uniform  temperature,  infinite  conductivity  model 
4b.  One-dimensional,  radial  temperature  distribution 

The  droplet  mass  is  continuously  updated  from: 

dm  _ 

~^  =  Spmv  +  mB  (5.7) 

Solving  the  energy  equation  for  the  droplet  is  a  boundary  condition  problem  and,  for 
models  4a  and  4b,  the  heat  flux  at  the  drop  surface  is  computed  first: 

Qs=~k% *'  +mv^vtt)  =  h(T„  -  Ts)  +  mvHv(Ts)  (5.8) 

Once  q.  is  known,  for  the  infinite  conductivity  model,  the  droplet  temperature  is  found  by  solving 
the  energy  equation  for  the  entire  drop: 


->■ 


(dTA 

i! 

pa 

k* 

(dPA 

k  dt  J 

Pv 

v  dt  j 

m„ 


(5.9) 


After  the  drop  rate  of  mass  and  temperature  change  have  been  updated,  the  rate  of  radius  change 
is  computed  as: 
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(5.10) 


dRp  =  \\  l  dmp  RpdPpdTp 
dt  pp  [4nR2p  dt  3  dT  dt 


with  the  implied  assumption  that  the  liquid  density  is  a  function  of  temperature  alone. 

To  obtain  the  radial  temperature  distribution,  the  one-dimensional  heat  conduction  equation 
is  written  in  spherical  non-dimensional  coordinates  and  solved  inside  the  drop: 


dP  l 


_ I _ & _ 


dRp 

dt 


dz 


dHD 


(S.ll) 


with  the  boundary  conditions: 


dT„ 


dz 


=  0,  q\z=l=Rpqs 


k=0 


(5.12) 


5.2.5  Liquid  Properties 

Earlier  studies  have  demonstrated  the  critical  importance  of  using  variable  properties  when 
computing  the  mass  and  heat  transfer  to  the  drop.  The  liquid  properties  are  presently  computed 
using  two  different  sets  of  approximate  formulae,  depending  upon  whether  or  not  the  liquid  is  an 
agent  simulant  or  standard  liquid.  When  the  liquid  is  a  simulant,  the  properties  are  computed 
according  to  the  curve  fits  provided  in  VLSTRACK.  For  other  liquids,  approximate  formulae  were 
obtained  from  Reid  [25],  and  these  are  also  listed  in  Appendix  B. 

5.2.6  Critical  Weber  Number  Criterion  Breakup  Models 

In  these  simple  models,  the  drop  diameter  is  continuously  reduced  to  prevent  the  drop 
Weber  number  from  exceeding  the  critical  Weber  number.  The  critical  Weber  number  is  computed 
as: 

We  =  Wec  + 12.924 Oh1  6  (5.13) 

where  the  Ohnesorge  number.  Oh  =  -r-L-  ,  provides  corrections  for  liquid  viscosity.  Equation 

^jpLd<j 

(5.13)  is  solved  iteratively  for  the  drop  diameter  at  each  time  step  and  the  drop  diameter  is  updated 
from: 
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dnew  nun  {d0i^ ,  d\ye ) 


(5.14) 


The  mass  loss  due  to  breakup  is  not  accounted  for  in  this  model.  The  implicit  assumption  is  that 
“daughter  droplets”  which  break  off  from  the  “parent”  are  small  enough  to  be  vaporized 
instantaneously. 

5.2.7  Dynamic  Deformation  and  Breakup  Model 

While  different  breakup  modes  occur  depending  on  the  Weber  number,  all  modes  start  with 
the  droplet  flattening  and  deforming  into  what  can  be  modeled  as  an  oblate  spheroid.  Modeling  the 
deformation  dynamically  makes  it  possible  to  track  the  shape  of  the  drop  as  it  moves,  and  therefore, 
to  take  into  account  the  drop  shape  in  the  drag  calculation  as  well  as  to  obtain  a  dynamic  criterion 
for  secondary  breakup.  We  are  currently  using  the  models  of  Ibrahim,  et  al.  [26],  and  Liu  and  Reitz 
[27]. 


In  their  model,  Ibrahim  [26]  assume  that  the  interaction  of  aerodynamic  forces  with  viscous 
and  surface  tension  forces  causes  the  drop  to  deform  into  an  oblate  spheroid.  The  energy  balance 
of  a  drop  under  purely  extensional  flow  is  written  and  yields  an  equation  for  the  oscillation  of  the 
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center  of  mass  of  a  half  drop  as  a  spring/mass/damper  system.  In  the  final  equation,  the  excitation 
is  provided  by  the  aerodynamic  pressure  forces,  the  surface-tension  forces  act  as  a  spring,  and  the 
viscous  forces  dissipate  the  energy.  In  the  present  formulation,  and  contrary  to  Ibrahim's  original 
model,  the  exact  expression  is  used  for  the  variation  of  surface  area  with  length  of  major  axis  after 
it  was  found  that  the  approximate  expression  proposed  by  Ibrahim  et  al.  yields  erroneous  values 
for  small  deformations.  The  procedure  is  as  follows. 

1 .  Solve  for  deformation  parameter,  y ,  (We  based  on  radius) 


dt  We  p  AR  da  p  y  dt  8  p.  R. 


(5.15) 


The  variation  of  surface  area  with  length  of  major  axis  is  given  by: 


dS, 


da 


p  _ 


=  Am  +  K 


Rb 


a 


.  .l  +  £,,  4  6 R  . 

]n(— ^)  + 


6R° 


1  —  £  '  e  2 aV  eV( l-ez) 


(5.16) 


2.  Calculate  oblate  spheroid  major  axis,  a 
3a.  Compute  new  frontal  area:  Af  =  II a2 
3b.  Modify  Cd  such  that 


Q>  =  Cd,  sphere)  1  +  2.632  min 


1, 


(5.17) 


Drop  breakup  is  assumed  to  occur  when  the  non-dimensional  deformation  parameter  equals  two, 
which  implies  that  the  length  of  the  spheroid’s  major  axis  is  twice  the  original  drop  radius.  Then, 
the  drop  breakup  is  computed  as  a  rate  process  using  correlations  for  breakup  times  [16]  and  for 
secondary  breakup  diameter  distribution  [27].  This  is  another  very  important  difference  with  the 
previous  simple  We  criterion  model,  where  the  mass  removed  by  breakup  is  simply  ignored.  Here, 
the  full  accounting  of  the  droplets  that  breakup  is  available. 

5.3  DEMONSTRATIVE  CALCULATIONS 

Four  series  of  results  that  highlight  the  advanced  features  of  SDROP  are  presented  next. 
For  all  the  results,  the  same  initial  conditions  were  used.  A  1  cm  diameter  VX-simulant  drop  is 
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released  at  90  km  above  sea  level  with  a  1,000  m/s  velocity  towards  the  center  of  earth  having  an 
initial  uniform  temperature  of  300  K.  The  four  different  cases  are:  1 -Uniform  T,  no  deformation,  2- 
T(r),  no  deformation,  3-Uniform  T,  simple  We  criterion  model,  and  4-Uniform  T,  dynamic 
deformation  model. 

Figure  5. 1  shows  the  drop  altitude  versus  time  for  the  four  cases  studied.  The  We  number 
criterion  results  are  represented  by  the  lighter  trace  which  shows  that  all  particles  are  broken-up 
(and  hence  are  removed  from  the  analysis)  within  50  s.  While  no  difference  is  noticeable  between 
the  results  without  deformation  (uniform  T  vs.  T(r)),  when  drop  deformation  is  taken  into  account, 
it  takes  approximately  twice  as  long  for  the  drop  to  reach  the  ground.  This  behavior  is  explained 
by  the  differences  in  drop  drag,  which  are  discussed  later  in  this  section.  The  time  evolution  of  the 
drop  diameter  is  plotted  in  Figure  5.2,  which  shows  that  the  particle  size  remains  almost  constant 
for  the  first  12  s  of  flight  and  then,  for  the  We  criterion  only,  decreases  abruptly. 

This  fundamental  difference  is  explained  by  displaying  the  We  number  versus  time  (Figure 
5.3).  In  the  We  criterion  calculations.  We  =  12  at  =  12  s,  and  remains  constant  thereafter  as  the 
drop  diameter  is  continuously  reduced.  As  the  drop  diameter  is  decreased,  the  heat  transfer  to  the 
drop  and  the  drop  temperature  increase,  until  the  drop  flash  vaporizes  and  the  calculation  is 
stopped.  In  the  deformation  results,  as  We  increases  (Fig.  5.3),  the  deformation  parameter 
increases  until  We  =  38  and  y  =  1.4,  at  which  time  both  We  and  y  start  decreasing. 

To  evaluate  the  influence  of  the  temperature  distribution  on  the  outcome  of  the  calculations, 
the  results  for  cases  1-3  are  displayed  in  Figure  5.4.  Even  though  the  absolute  change  in  drop 
diameter  is  minor  in  all  cases,  resolving  the  radial  temperature  distribution  is  seen  to  yield  a  sharp 
increase  in  vaporization.  The  increased  evaporation  is  directly  related  to  the  higher  surface 
temperature  of  the  drop,  as  shown  in  Figure  5.5.  As  mentioned  earlier,  the  rate  of  surface 
temperature  rise  is  greatest  for  the  We  criterion  model  due  to  the  shrinking  size  of  the  droplet, 
which  promotes  flash  vaporization  by  yielding  hotter  drops  at  higher  altitudes.  The  lowest  surface 
temperature  is  found  for  the  deformation  case,  where  the  drop  is  slowed  by  the  excess  drag.  The 
temperature  distribution  inside  the  droplet  for  case  2  is  shown  at  different  times  in  Figure  5.6.  The 
surface  temperature  first  increases  (curve  B)  and  then  decreases  (curves  C-I)  after  the  drop  slows 
down.  As  the  drop  reaches  the  ground,  the  air  temperature  increases  and  the  drop  temperature 
increases  again. 

The  significant  differences  in  Co  and  drop  velocity  for  cases  1,2  and  3  are  illustrated  in 
Figures  5.7  and  5.8.  From  Figure  5. 1,  the  first  20  km  of  flight  are  completed  in  about  50  s  and  the 
sharp  increase  in  drop  drag  (Fig.  5.7)  is  directly  correlated  to  the  drop  deformation  (Fig.  5.3).  At 
comparable  altitudes,  the  deformed  drop  velocity  is  considerably  lower  than  without  deformation 
(Fig.  5.8),  which  explains  the  large  discrepancy  in  total  flight  time  obtained  in  Figure  5.1.  An 
important  consequence  of  this  behavior,  which  is  not  illustrated  in  these  calculations,  is  that  the 
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increased  Cd  and  flight  times  would  yield  landing  locations  that  are  very  different  from  those 
without  deformation.  The  drop  shape  history  is  provided  in  Figure  5.9.  In  the  first  30  s  of  flight 
the  drop  flattens,  then  relaxes  as  the  velocity  and  the  aerodynamic  pressure  decrease.  In  the  lower 
and  denser  reaches  of  the  atmosphere,  the  deformation  increases  again  to  reach  the  final  shape  at  Z 
=  0  km. 

The  example  presented  indicates  that  with  use  of  a  more  advanced  breakup  model, 
survivability  of  droplets  is  greatly  enhanced,  whereas  the  simple  critical  Weber  number  criterion 
model  indicates  that  none  survive.  Emphasis  must  be  placed  upon  use  of  advanced  models  for  this 
problem  as  demonstrated  by  this  example. 
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Fig.  5.2.  Drop  diameter  versus  time. 


-46- 


40 


1.5 


Fig.  5.3.  We  number  (solid)  and  deformation  parameter  (dash)  versus  time. 
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Fig.  5.7.  Drop  Cd  versus  altitude. 
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Fig.  5.9.  Cross  sections  of  the  drop  at  different  points  on  the  trajectory. 


6.0  ADVANCED  GENERALIZED  MULTIPHASE  FORMULATION  ON 

UNSTRUCTURED  GRIDS 


6.1  OVERVIEW 

The  CRAFT  code  has  direct  applicability  to  post-hit  problems  involving  the  survivability  or 
neutralization  of  droplet  (or  particulate)  clouds,  but  limitations  in  addressing  aspects  of  the  primary 
breakup  problem  from  a  first  principles  perspective.  The  breakup  problem  requires  numeric  which 
can:  (1)  Concurrently  analyze  nearly  incompressible  liquids  and  highly  compressible  gases 
(maintaining  acoustic  accuracy,  since  it  is  the  instability  at  the  gas/liquid  that  leads  to  breakup);  and, 
(2)  adapt  the  grid  dynamically  to  the  changing  structure  of  the  gas/liquid  interface,  maintaining 
sufficient  resolution  to  permit  a  first-principles  treatment  of  breakup  for  resolvable  droplet  sizes. 
Presently,  such  numerics  are  not  available,  with  research  required  in  the  areas  of  multi-phase 
preconditioning  and  multi-phase  unstructured  numerics.  Another  feature  required  is  that  of 
cavitation  to  address  breakup  related  to  the  rapid  escape  of  gas  bubbles  within  a  liquid  glob  at 
higher  altitudes.  We  discuss  below  some  background  and  the  path  taken  to  evolve  these  capabilities 
derived  from  related  activities. 

6.2  BACKGROUND 

Over  the  years,  the  diverse  requirements  in  analyzing  fluid  flow  regimes  with  varying 
features  has  led  to  independent  algorithmic  development  and  modeling  initiatives  for  each  class  of 
problems.  In  the  incompressible  limit,  pressure  based  methods  were  developed  to  solve  for  the 
momentum  equations  of  the  Navier-Stokes  set  and  a  Poisson  solver  was  used  to  determine  the 
pressure  field.  In  the  compressible  limit,  the  Navier-Stokes  system  could  be  modeled  as  a  complete 
hyperbolic  system  and  all  variables  could  be  solved  in  a  unified  manner  using  density-based 
methods  (such  as  in  CRAFT).  However,  in  our  problem,  multiple  flow  regimes  and  fluid  phases 
need  to  be  analyzed  concurrendy  which  creates  an  inherent  need  to  solve  for  the  different  regimes 
of  flows  in  a  consistent  manner  within  the  same  general  framework. 

Chorin's  artificial  compressibility  method  [28]  initiated  the  unification  process  by  casting 
the  incompressible  equations  in  a  hyperbolic  manner  through  the  use  of  a  pressure/pseudo-time 
derivative  in  the  continuity  equation.  As  a  consequence,  the  form  of  the  incompressible  equation  set 
is  hyperbolic  and  it  reduces  to  the  incompressible  form  in  the  steady  state.  Calculation  of  unsteady 
incompressible  problems  are  carried  out  by  using  a  dual-time  stepping  procedure  [29],  where  the 
time  derivative  term  is  included  as  a  source  term  and  at  each  time  step,  the  solution  is  driven  to 
convergence  based  on  pseudo-time  stepping.  Preconditioning  techniques  [30,31]  have  been 
developed  to  provide  a  transformation  that  enables  the  solution  of  all  flow  regimes  from  the 
incompressible  limit  to  the  highly  compressible  limits.  The  preconditioned  system  utilizes  Chorin's 
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concept  of  artificial  compressibility  and  takes  it  a  step  further  by  redefining  certain  terms  of  the 
transformation  matrix  with  the  help  of  asymptotic  analysis  and  arguments  based  on  order  of 
magnitude.  This  greatly  alleviates  convergence  problems  by  rescaling  the  eigenvalues  especially  in 
the  incompressible  limit  where  the  equations  become  stiff. 

Preconditioned  methods  have  been  shown  to  be  robust  for  single-phase  problems  ranging 
from  incompressible  type  flows  to  highly  compressible  flows.  However,  these  formulations  seldom 
cover  two-phase  flow  regimes,  where  a  highly  compressible  flow  regime  coexists  with  an 
incompressible  regime.  The  physical  phenomena  observed  in  these  multiphase  flows  are  dictated 
by  the  interaction  between  the  two  phases.  Small  changes  in  interfacial  properties  can  manifest 
themselves  as  significant  alterations  of  the  macroscopic  properties  of  individual  phases.  Therefore, 
accurate  modeling  of  the  interfacial  dynamics  of  such  systems  becomes  critical.  A  class  of  methods 
called  interface  fitting  methods  has  been  used  explicitly  to  track  and  fit  the  gas/liquid  interface  with 
an  emphasis  on  the  treatment  of  cavitating  flows.  However,  these  methods  [32-37]  utilize  potential 
flow  and  closed  cavity  approximations  that  are  not  appropriate  for  the  problem  of  interest  where  the 
interface  deforms  substantially. 

Interface  capturing  methods  as  used  in  CRAFT  are  attractive,  in  part,  because  the  gas/liquid 
interface  is  an  outcome  of  the  solution  procedure.  As  a  consequence,  grid  resolution  and  the 
numerical  algorithm  strongly  affect  solution  quality  at  the  gas-liquid  interface.  This  places  the 
burden  on  the  closure  model  and  the  numerics  to  accurately  reflect  the  thermodynamics  associated 
with  the  two-phase  system.  Furthermore,  numerical  diffusion  can  exhibit  a  smearing  effect  since 
the  interface  is  an  artifact  of  the  solution.  This  problem  could  potentially  be  resolved  by  increased 
grid  resolution  locally  close  to  the  interface.  However,  the  interface  is  not  known  apriori  impeding 
efforts  to  provide  the  requisite  amount  of  grid  resolution  in  the  appropriate  region  of  the 
computational  domain.  Recently,  in  work  performed  by  CRAFT  Tech  carried  out  by  Cavallo  and 
Baker  [38],  a  novel  method  of  grid  adaption  for  unstructured  mixed  element  meshes  has  been 
developed.  In  this  work,  we  have  evolved  the  capability  to  refine  the  grid  based  on  gradients  of  local 
void  fraction  by  a  Delaunay  cavity  reconstruction  based  adaption  procedure  [38]. 

In  using  unstructured  numerics  with  dynamic  adaption,  we  bridge  the  gap  between  interface 
tracking  methods  and  interface  capturing  methods,  as  the  interface  is  part  of  our  solution  but 
successive  adaption  refines  it  as  the  simulation  progresses.  This  procedure  leads  to  savings  in 
computational  resources,  since  the  interface  is  first  identified  and  then  the  grid  in  the  vicinity  of  the 
interface  is  refined.  In  work  described  below,  we  present  a  generalized  mixture  formulation  model 
that  extends  the  framework  of  our  CRUNCH  compressible  unstructured  code  [39,40]  to  include 
varied  elements  required  for  the  liquid  breakup  problem.  We  describe  this  in  progressive  steps. 
The  first  subsection  contains  the  derivation  for  the  multi-phase  mixture  formulation.  We  initially 
consider  two  constant-density  fluids,  with  constant  but  different  acoustic  speeds.  This  is  important 
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to  demonstrate  the  accuracy  of  the  code  and  the  formulation  in  the  incompressible  limit.  It  is  also 
important  because  we  show  that  the  acoustics  of  the  multiphase  system  is  consistent  with  the 
thermodynamics  of  the  system.  Here  the  local  speed  of  sound  in  the  two-phase  mixture  is  a 
function  of  the  local  void  fraction  and  closely  mimics  the  two-phase  acoustic  speed  relationship 
from  classical  analytic  theory.  This  is  followed  up  by  two  cases  for  validation  purposes.  In  the 
following  section,  we  extend  the  mixture  formulation  to  include  compressibility  effects.  The  code  is 
now  operational  in  all  flow  regimes,  and,  for  both  the  incompressible  and  compressible  phases,  with 
an  equation  of  state  being  provided  for  the  case  of  each  phase. 

6.3  MULTIPHASE  EQUATION  SYSTEM 
6.3.1  Original  Multiphase  Equations 

The  multiphase  equation  system  is  written  in  vector  form  as: 


» 


dQ  dE  dF  dG  =  v 
dt  dx  +  dy  +  dz  x'y' 


+  S 


(6.1) 


Here  Q  is  the  vector  of  dependent  variables,  E,  F  and  G  are  the  flux  vectors,  S  the  source  terms,  and 
^xyz  represents  the  viscous  fluxes.  The  details  of  the  source  term  S  for  a  cavitation  process  are 
discussed  later  in  the  section.  The  vectors  Q  and  E  are  given  below  and  should  be  compared  with 
those  of  the  multi-phase  CRAFT  code  as  given  by  Eqn.  2.2.  CRAFT  contains  an  energy  (em) 
equation  and  supplemental  equations  for  gaseous  species/combustion  products  not  in  the  initial 
incompressible  CRUNCH  formulation.  CRUNCH  shows  equations  for  turbulence  variables  k  and 
£  which  are  also  in  CRAFT  but  not  which  shown  (for  clarity)  in  Eqn.  2.2. 
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The  vectors  F  and  G  can  be  written  similarly.  pm  is  the  mixture  density,  and  <j>g  is  the  volume 
fraction  or  porosity  for  the  gas  phase.  Note  that  in  Eqn.  (6.2)  the  energy  equation  is  dropped  since 
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each  phase  is  assumed  to  be  nearly  incompressible  thereby  decoupling  the  pressure  work  term. 
The  mixture  density  and  gas  porosity  are  related  by  the  following  relations: 


Pm  ~  PgPg  +  PlPl 

(6.3) 

^  =  Pg+^L 

(6.4) 

where  pg,  pL  are  the  physical  material  densities  of  the  gas  and  liquid  phase  respectively. 


6.3.2  Incompressible  Modification 

To  modify  the  system  in  Eqn.  (6.1)  to  a  well-conditioned  form  in  the  incompressible  regime 
requires  a  two-step  process;  an  acoustically  accurate  two-phase  form  of  Eqn.  (6.1)  is  first  derived, 
followed  by  a  second  step  of  time-scaling  or  preconditioning  to  obtain  a  well-conditioned  system. 
We  begin  by  defining  the  acoustic,  differential  form  of  the  mixture  density  pm  as  follows: 


dpm=(pg-pL)d<t>g+-\dP 

c<p 

„  2  „  2  T  „  2 


(6.5) 


is  a  variable  defined  for  convenience  and  is  not  the  acoustic  speed,  Cm,  in  the  mixture  which  will 


be  defined  later.  Here  c  is  the  isothermal  speed  of  sound 


dP 


l  dP, 


in  the  pure  gas  phase,  and  cL  is  the 


corresponding  isothermal  speed  of  sound  in  the  liquid  phase,  which  is  a  finite-value.  The  original 
system  of  equations  (6.1)  is  transformed  by  replacing  the  dependent  vector  array  Q  by  the  vector  of 
primitive  variables  Qv.  There  are  many  inherent  advantages  in  solving  the  system  numerically  for 
the  primitive  variable  set  Qv.  The  pressure  field  is  now  directly  integrated  as  part  of  the  solution 
vector.  This  is  especially  beneficial  when  solving  for  incompressible  type  flows.  The  recasting  of 
the  solution  vector  also  makes  the  system  more  amenable  to  implementation  of  convergence 
acceleration  techniques  such  as  preconditioning,  since  the  preconditioning  matrix  is  a  modification 
of  the  transformation  matrix,  Gamma  [T].  The  viscous  fluxes  are  also  more  readily  computed  since 
the  derivatives  of  the  primitive  variable  set  are  involved  in  their  calculation.  Using  Eqn.  (6.5),  Eqn. 
(6.1)  may  be  rewritten  as 


dQ,  dE  dF  dG 
dt  dx  dy  dz 


=  V 


+  S 


(6.6) 
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The  numerical  characteristics  of  Eqn.  (6.7)  are  studied  by  obtaining  the  eigenvalues  of  the  matrix, 
which  are  derived  to  be: 


p— 1 

f  <3E  ) 
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teJ. 

A  =  (u  +  cm,  u-cm,  u,u,u,u,u)  (6.8) 

where  cm  turns  out  to  be  well-known,  harmonic  expression  for  the  speed  of  sound  in  a  two-phase 
mixture  and  is  given  as: 
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The  behavior  of  the  two-phase  speed  of  sound  is  plotted  in  Fig.  (6.1)  as  a  function  of  the 
gas  porosity;  at  either  limit  the  pure  single-phase  acoustic  speed  is  recovered.  However,  away  from 
the  single-phase  limits,  the  acoustic  speed  rapidly  drops  below  either  limit  value  or  remains  at  the 
low-level  in  most  of  the  mixture  regime. 

We  reemphasize  a  critical  observation  at  this  point:  the  Equation  system  6. 6-6.7  is 
completely  defined  and  does  not  require  ad-hoc  closure  models  for  the  variation  of  mixture  density 
with  pressure.  In  that  respect  alone  this  represents  a  significant  advancement  over  most  other 
multiphase  models  in  the  literature.  The  acoustic  speeds  for  individual  phases  are  well-defined 
physical  quantities,  which  may  be  specified,  and  so  is  the  case  with  physical  material  densities  (pg. 
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pL)  for  each  individual  phase.  For  low  pressure  incompressible  regimes,  the  material  densities  may 
be  assumed  to  be  constant  without  significant  error  in  the  solutions.  However,  in  its  most  general 
form  the  material  densities  for  each  phase  may  be  obtained  from  the  pressure  using  their  respective 
physical  equations  of  state  (e.g.  ideal  gas  law  for  gases,  etc)  if  that  is  so  desired,  and  will  be  taken 
up  later  in  this  section. 


6.3.3  Preconditioning  of  Equation  System 

Preconditioning  is  now  applied  to  the  system  in  an  effort  to  rescale  the  eigenvalues  of  the 
system  so  that  the  acoustic  speeds  are  of  the  same  order  of  magnitude  as  the  local  convective 
velocities.  This  is  achieved  by  replacing  T  in  Eqn.  (6.6)  by  T . 
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Here,  the  parameter  (3  has  been  introduced  to  precondition  the  eigenvalues.  The  modified 
eigenvalues  of  the  preconditioned  system  are  given  as: 
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where 
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Eqn.  (6.12)  indicates  that  by  setting  /3  =  ^  2  J  where  up  =  max  (u,  0.01  cm)  the  pseudo-acoustic 

speed  is  of  the  order  of  u  at  all  mixture  composition  values.  Therefore,  a  rigorous  definition  of  the 
physical  two-phase  acoustic  speed  cm  (as  has  been  done  here)  is  critical  to  obtaining  noise-free 
propagation  of  pressure  waves  across  interfaces  where  the  density  and  acoustic  speed  are  varying 
rapidly. 


6.3.4  Cavitation  Extensions 

For  nearly  incompressible  gas/liquid  systems,  cavitation  problems  can  serve  to  validate  the 
numerics.  A  simplified  finite-rate,  rate  constant  as  applied  to  cavitation  problems  can  be  specified 
as  follows: 


wg=KepL(l>L  +  Kvpg<l>g 


(6.13) 


where  the  constants  Kx  provides  the  rate  constant  for  vapor  being  generated  from  liquid  in  a  region 
where  the  local  pressure  is  less  than  the  vapor  pressure.  Conversely,  Kx  provides  the  reconversion 
of  vapor  back  to  liquid  in  regions  where  the  pressure  exceeds  the  vapor  pressure.  While  we  discuss 
here  the  development  of  the  source  terms  for  cavitation  type  problems,  appropriate  finite  rate  source 
terms  can  be  considered  for  most  problems  undergoing  phase  change  due  to  vaporization,  etc. 
Here,  the  rate  constants  are  specified  using  the  form  given  below: 


K(  = 


0 

1 


P<PV 

h£z£H  p>pv 


CT «  \p~Q~ 


c  =  10"3,  x  =  Cav.  No.  =  Pv 
<2L  1  ~ 1 


(6.14) 


We  note  that  the  finite-rate  constant  has  not  been  correlated  with  experimental  data.  For 
steady  attached  cavitation,  this  simplified  form  may  be  adequate  since  the  cavitation  time  scales  do 
not  interact  with  the  fluid  time  scales  eventually  providing  the  equilibrium  conditions  at  steady  state. 
For  unsteady  cavitation,  however,  the  details  of  how  the  non-equilibrium  source  term  is  specified 


-57- 


could  be  crucial  since  it  may  couple  with  transient  pressure  waves.  The  development  of  a  more 
rigorous  non-equilibrium  source  term  model  is  a  topic  for  ongoing  inquiry. 


6.3.5  Numerical  Flux  Construction 

The  numerical  flux  at  each  cell  interface  is  obtained  by  solving  the  Riemann  problem 
locally.  Roe’s  flux  difference  method  is  utilized  to  construct  the  approximate  solution  to  the 
Riemann  problem.  The  flux  at  the  interface  is  constructed  by  evaluating  the  change  in  flux 
associated  with  each  wave  component  characterized  by  the  eigenvalues  of  the  Jacobian  matrix. 

[r-'A„],  where 
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.  The  flux  at  the  interface  is  given  by 

f,=j{fl+fs+af:+*f;)  (6.15) 

The  eigensystem  analysis  pertaining  to  [r'X]  reveals  the  eigenvalues  of  the  system  as 

(u  +  cm,  u-cm,  u,u,u,u,u)  as  denoted  by  Eqn.  (6.8).  The  right  and  left  eigenvectors  associated 
with  [r-'\]  are  given  as 


R  = 


0 


0 


psPm 

P*Cp. 

4 

K 

PSo. 

Pm^pm 

1 

4 

PSPm 

P"Ppm 

A 

A 

0 


0 


0  0  0  0  0 

mx  nx  ooo 

my  ny  0  0  0 

ml  nz  0  0  0 

0  0  10  0 

0  0  0  1  0 

0  0  0  0  1 


(6.16) 


-58- 


L  = 


K 

2  KPm~pm 

2  lyPSPm 

2  KPmCpm 

~\lxPrnCPm 

~2^yP m^Pm 

~~zhPmCPm 

0 

my 

mz 

0 

ny 

nz 

A 

PSl 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

in  Eqn.  (6. 15)  can  be  written  as 
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where  a  represents  the  matrix  of  eigenvalues . 

R,  L,  and  a  are  computed  at  the  interface  using  their  Roe-averaged  values.  The  dependent 
variables  that  comprise  these  matrices  are  computed  as 
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The  void  fraction  at  the  interface  is  determined  by  taking  an  arithmetic  mean  of  the  void  fractions  of 
the  left  and  right  states. 
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6.3.6  Inclusion  of  Compressible  Flow  Effects 

So  far,  we  have  derived  the  multi-phase  formulation  for  two  fluids  with  constant  densities. 
In  this  system  the  energy  equation  is  decoupled  from  the  momentum  equations.  Next,  we  derive  the 
multiphase  system  formulation  including  compressibility  effects.  In  this  derivation,  we  assume  an 
equilibrium  formulation  i.e.  we  will  work  with  one  temperature  for  both  fluids  in  a  gas/liquid 
computational  cell.  Each  fluid  can  have  a  user-defined  equation  of  state.  In  the  case  of  a 
liquid/incompressible  phase,  a  stiff  equation  of  state  can  be  specified.  The  framework  is  general 
and  can  accommodate  any  user  specified  equation  of  state. 

As  before,  the  governing  Navier-Stokes  equations  along  with  a  species  transport  equation 
can  be  expressed  in  the  form 


r%k+te  dF 

dt  dx  dy  dz 


=  Vxy. 


+  S 


where  now  the  primitive  variable  set  includes  temperature 


(6.24) 


T 

Qv  =  [p,u,v,w,(j)g,T^  (6.25) 

Again,  as  before,  the  functional  relationship  that  relates  the  mixture  density  with  the  local  gas 
porosity  remains  the  same 


Pm=Pg<t>g+PL<l>L  (6-26) 

1  =  <f>g+<PL  (6-27) 

where  pg,  pL  are  the  physical  material  densities  of  the  gas  and  liquid  phase  respectively.  The 
pressure  at  any  point  in  the  mixture  can  be  related  to  the  local  temperature  and  liquid/gas  densities 
by  the  user  defined  equation  of  states  for  the  liquid  and  gas  respectively.  As  a  consequence,  the  T 
matrix  indicative  of  the  conservative  variables  derivative  with  respect  to  the  primitive  variables  can  be 
represented  as: 
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v  j 

where  the  specific  heat  for  the  mixture  can  be  based  on  the  specific  heats  for  the  gas  and  the  liquid 
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Similarly,  the  enthalpy  for  the  mixture  can  be  defined  as  a  linear  combination  of  the  enthalpy  for  the 
liquid  and  the  gas 
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and  for  both  the  liquid  and  the  gas 
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The  numerical  characteristics  of  the  system  are  studied  by  obtaining  the  eigenvalues  of  the  matrix, 
which  are  derived  to  be: 
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where  c  turns  out  to  be  the  well-known,  harmonic 


expression  for  the  speed  of  sound  in  a  two-phase  incompressible  mixture  and  is  given  as  before  as 
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It  should  be  noted  that  in  the  limit  of  a  mixture  of  two  purely  incompressible  liquids,  the  energy 
equation  decouples  and  the  acoustic  speed  in  the  mixture  reverts  to  the  above  equation.  The  left  and 
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6.4  VALIDATION  STUDIES 

6.4.1  Single  Phase  Validation  -  Incompressible  Flow  in  a  90°  Square  Duct 

A  suitable  test  case  is  the  computation  of  laminar  flow  through  a  90  degrees  curved  square 
duct.  This  test  case  has  been  well  studied,  experimentally  verified  and  documented  by  Humphrey  et 
al.  [41]  and  has  historically  served  as  a  notable  benchmark  case  as  being  representative  of  some 
small  measure  of  flow  complexity.  Most  notably,  the  case  exhibits  secondary  flows  that  are  set  up 
in  the  bend  region  due  to  vorticity  being  generated  at  the  entry  to  the  bend,  thereby  forcing  the 
maximum  streamwise  velocity  in  the  plane  of  symmetry  towards  the  outer  wall.  The  Reynolds 
number  for  the  case  is  790  (based  on  the  hydraulic  diameter).  The  mesh  for  this  case  consisted 
entirely  of  hexahedral  elements  and  was  clustered  towards  the  walls  of  the  duct.  The  grid  generated 
for  the  problem  extended  10  diameters  upstream  and  downstream  of  the  bend. 

Figure  6.2  shows  the  pressure  distribution  on  the  plane  of  symmetry.  The  figure  clearly 
shows  the  higher  pressure  region  setting  up  close  to  the  outer  wall  near  the  bend  as  a  consequence 
of  the  turning  of  the  flow.  Figure  6.3  shows  the  streamwise  velocity  comparison  between  the 
experimental  observations  and  the  simulation  at  the  exit  of  the  bend.  The  CFD  calculation  captures 
the  velocity  variation  close  to  the  inner  and  outer  walls  and  is  in  fairly  good  agreement  with  the 
experimental  measurements. 

6.4.2  Multiphase  Validation  -  NACA  66  Hydrofoil  with  Cavitation 

The  performance  of  marine  propellers  is  severely  affected  when  the  blades  have  to  operate 
under  off-design  conditions,  hi  particular,  these  blades  are  subject  to  unsteady  pressure  oscillation 
in  turbulent  wakes.  Under  these  circumstances,  the  local  pressure  close  to  the  surface  of  the  blade 
has  a  tendency  to  go  below  the  vapor  pressure  giving  rise  to  the  formation  of  cavities/vapor  pockets. 
This  significantly  increases  the  drag  on  the  propeller  and  places  a  limitation  on  the  thrust  generated 
by  the  propeller.  As  a  result,  simulations  of  cavitation  over  hydrofoils  become  important  in  aiding 
the  understanding  of  the  inception  mechanisms  and  the  dynamics  between  bubbles  and  the  pressure 
fluctuations  on  the  blade  surface. 

The  hydrofoil  being  utilized  in  this  test  case  is  a  modified  NACA  66  hydrofoil,  that  is 
commonly  utilized  in  propeller  design.  Shen  and  Dimotakis  Ref.  [42]  have  documented 
experimental  observations  related  to  cavitation  on  the  NACA  66  hydrofoil.  The  grid  constructed  to 
study  this  case  consists  of  a  combination  of  194,751  tetrahedral  and  prismatic  elements.  The  grid 
was  tightly  clustered  close  to  the  surface  of  the  hydrofoil  especially  in  the  upper  surface  region  of 
the  hydrofoil.  The  Reynolds  number  for  this  case  was  2  x  106.  A  density  and  viscosity  ratio  of 
1:100  was  maintained  between  the  two  phases.  A  flow  at  an  angle  of  attack  of  4  degrees  was 
modeled  at  cavitation  numbers  of  0.89  and  0.91. 
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Pressure  contours  for  the  hydrofoil  at  a  cavitation  number  of  0.84  are  plotted  in  Figure  6.4. 
This  initial  calculation  was  performed  with  a  wall-function  procedure.  Note  that  the  pressure 
contours  cluster  around  the  cavitation  boundary  and  resemble  a  weak  shock  front.  Since  the 
acoustic  speed  drops  significantly  in  the  interface  region,  the  flow  transitions  from  subsonic  to 
supersonic  in  the  interface  and  back  to  subsonic  inside  the  cavity  giving  rise  to  the  clustering  in  the 
pressure  contours.  The  gas  void  fraction  contours  on  the  hydrofoil  surface  are  shown  in  Figure  6.5 
for  a  cavitation  number  of  0.9.  The  cavitation  is  enclosed  within  a  very  narrow  region  up  to  30 
percent  of  the  chord  highlighting  the  need  for  a  grid  adaption  procedure  which  will  locally  refine 
cells  around  the  interface.  The  surface  pressure  distribution  is  compared  with  experimental  data  in 
Figure  6.6  for  this  case.  We  observe  that  the  predicted  cavity  is  smaller  than  the  experimental 
measurement.  Also  the  pressure  recovery  at  the  mid-chord  is  not  very  accurate.  These  deficiencies 
were  attributed  to  the  wall  function  procedure  which  is  well  known  to  underpredict  separated  flow 
regions. 

In  order  to  determine  the  relative  performance  of  a  wall  function  procedure  versus  a  near¬ 
wall  turbulence  model,  we  computed  the  flowfield  at  a  cavitation  number  of  0.84  using  both  models. 
The  void  fraction  contours  for  this  case  are  shown  in  Figure  6.7.  The  characteristics  of  the  cavity 
are  similar  to  the  previous  case  except  that  the  cavity  extends  to  about  50  percent  of  the  chord  due  to 
the  lower  cavitation  number.  The  surface  pressure  profile  is  plotted  in  Figure  6.8  with  the 
calculations  done  with  both  a  wall  function  as  well  as  a  near  wall  turbulence  model  of  So  and  Zhang 
[43].  The  wall  function  procedure  as  before  overpredicts  the  aft  recovery  pressure  significantly. 
However,  with  a  near- wall  turbulence  model  good  comparisons  are  obtained  for  the  pressure  behind 
the  cavity.  The  length  of  the  cavity  does  not  show  much  sensitivity  with  both  models  having 
approximately  the  same  value  and  that  is  close  to  the  experimental  measurement. 

6.4.3  Single  Phase  Compressible  Validation  -  Supersonic  Ramp 

The  multiphase  formulation  has  been  successfully  validated  at  the  incompressible  limit.  The 
next  step  is  to  use  the  multiphase  formulation  for  validation  at  the  compressible  limit  for  a  single 
phase.  The  case  considered  for  this  is  supersonic  flow  over  a  ramp  at  an  approach  Mach  number  of 
3.0.  This  is  a  standard  test  case  and  we  compare  the  solution  of  the  multiphase  CRUNCH  code 
with  that  obtained  from  the  standalone  compressible  CRUNCH  code  [39,40].  The  comparisons  are 
shown  in  Figure  6.9,  where  pressure  contours  are  plotted  from  both  solutions.  The  multiphase 
solution  compares  favorably  with  the  compressible  code  solution.  The  series  of  compressions  and 
expansions  are  captured  very  accurately.  The  compression/expansion  wave  locations  at  the  wall  are 
almost  identical.  The  strength  of  the  compressions  from  the  multiphase  solutions  also  compares 
very  well  with  the  compressible  solution  with  very  little  or  no  smearing. 
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6.4.4  Compressible  Multiphase  Two-Phase  -  Shock  Tube  Studies 

Two  series  of  shock  tube  cases  were  computed  with  a  pressure  ratio  of  2:1  with  the  high 
pressure  side  at  lOOMPa  and  the  low  pressure  side  at  50  MPa.  In  the  first  series,  the  liquid  starts 
out  on  the  low  pressure  side  and  the  gas  is  at  the  high  pressure  side.  In  the  next  series,  the  gas  is 
subjected  to  the  lower  pressure  and  the  liquid  starts  out  on  the  higher  pressure  side.  In  both  cases, 
we  have  captured  the  correct  shock  amplitude  and  the  propagation  of  the  shock  and  rarefaction 
waves.  Figure  6.10  shows  the  solutions  for  the  two  series  at  various  snapshots  in  time. 

6.4.5  Low  speed  Bulk  Liquid  Flyout  Deformation  in  Air 

This  test  case  consists  of  an  analysis  of  deformation  of  a  liquid  blob  without  breakup.  The 
blob  simulations  performed  in  Section  3  with  the  CRAFT  code  were  all  in  the  supersonic  regime. 
With  the  new  multiphase  formulation  based  on  the  primitive  variables,  we  can  now  look  at  blob 
deformations  in  the  low  Mach  number  regimes,  as  the  flow  physics  and  the  physics  of  deformation 
here  will  be  very  different  from  that  encountered  earlier  in  Section  3. 

With  this  capability,  we  can  now  cover  all  regimes  of  blob  deformation,  from  the  high 
supersonic  to  the  low  subsonic.  The  grid  used  in  this  simulation  is  shown  in  Figure  6.1 1.  The  grid 
is  colored  by  void  fraction,  with  the  square  in  the  middle  of  the  computational  domain  indicative  of 
the  initial  shape  of  the  blob.  The  grid  is  clustered  close  to  the  gas-liquid  interface.  As  the  blob 
deforms,  it  goes  through  two  distinct  evolutionary  stages.  The  first  stage  is  shown  in  Figure  6.12 
where  the  blob  starts  flattening  and  arcing  outwards  forcing  the  flow  around  it.  This  leads  to  the 
formation  of  recirculating  regions  in  the  shadow  of  the  interface  arc  at  the  front  end  of  the  blob. 
Flow  behind  the  blob  behaves  like  the  flow  in  the  wake  of  a  blunt  body.  The  second  evolutionary 
stage  in  the  deformation  of  the  blob  (Fig.  6.13)  shows  a  distinct  flattening  of  the  blob,  with  the  blob 
getting  elongated  to  almost  twice  its  initial  height.  The  recirculation  region  behind  the  arc  at  the 
head  of  the  blob  has  considerably  shrunk  in  size.  There  is  a  sizeable  compression  of  the  liquid 
mass  in  the  flow  direction,  due  to  the  lateral  elongation.  The  lateral  elongation  further  increases  the 
size  of  the  recirculation  zones  at  the  aft  end,  as  the  oncoming  flow  now  has  to  negotiate  a  larger  area 
to  get  around  the  blob.  The  distinct  difference  with  the  high  speed  blob  discussed  in  Section  3  is 
the  elongation  of  the  blob  in  the  low  speed  case  and  the  absence  of  fingerlike  projections  that  are 
characteristic  of  blob  deformation  at  high  speeds. 
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Fig.  6.2.  Pressure  distribution  along  plane  of  symmetry  in  the  90°  degree  curved  square  duct. 
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Normalized  Distance 

Fig.  6.3.  Comparison  of  the  streamwise  velocity  along  the  plane  of 
symmetry  with  experimental  data  of  Humphrey  et  al.  [41]. 
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Fig.  6.4.  A  representative  pressure  distribution  on  the  NACA  66  hydrofoil  at  4  degrees  angle  of 

attack  and  cavitation  number  of  0.84. 
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Cavitation  No.  =  0.91  Reynolds  No.  =  2,000,000  Angle  of  Attack  =  4  deg 


Fig.  6.5.  Cavitation  bubbles  indicated  by  void  fraction  contours  on  the  NACA  66  modified 

hydrofoil  at  cavitation  numbers  of  0.91. 
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Fraction  of  chord 

Fig.  6.6.  Cp  comparisons  with  experimental  data  of  Shen  and  Dimotakis  [42]  for  the  NACA  66 

hydrofoil  at  a  cavitation  number  of  0.91. 


Cavitation  No.  =  0.84  Reynolds  No.  =  2,000,000  Angle  of  Attack  =  4  deg 


Fig.  6.7.  Cavitation  bubbles  indicated  by  void  fraction  contours  on  the  NACA  66  modified 

hydrofoil  at  cavitation  numbers  of  0.84. 
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O  CFD  -  Wall  function 


Fig.  6.8.  Cp  comparisons  with  experimental  data  of  Shen  and  Dimotakis  [42]  for  the  NACA  66 

hydrofoil  at  cavitation  number  of  0.89. 
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Fig.  6.9.  Flow  over  supersonic  ramp. 
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Fig.  6.10.  Pressure  traces  at  various  snapshots  for  the  shock  tube  cases. 


Yellow  region  indicates  gas  and  the  blue  region  is  the  initial  liquid  blob. 
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zones  behind  the  gas-liquid  front  and  in  the  wake  region. 


blob  at  low  speeds  and  the  sizeable  recirculation  patterns  in  the  wake  region. 
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7.0  CONCLUDING  REMARKS 


A  number  of  advanced  capabilities  have  evolved  from  the  work  performed  on  this  Phase  II 
SBIR  program.  The  CRAFT  code  has  demonstrated  a  rudimentary  capability  to  analyze  post-hit 
liquid  blob  deformation  and  breakup  processes  but  is  limited  in  its  ability  to  highly  resolve  the 
gas/liquid  interface  and  to  concurrently  analyze  the  interactions  between  highly  compressible  gases 
and  nearly  incompressible  liquids.  Exploratory  work  with  a  new  gas/.liquid  variant  of  the 
CRUNCH  unstructured  code  has  exhibited  the  potential  to  remedy  these  limitations  and  progress  is 
being  made  with  this  new  version  of  CRUNCH  to  deal  with  problems  in  cavitation  and  related 
multi-phase  problems. 

The  CRAFT  code  can  perform  quite  well  in  post-hit  scenarios  where  the  initial  droplet  (or 
particulate)  distribution  can  be  defined  to  initialize  the  solution.  CRAFT  is  currently  being  used  to 
analyze  bunker  neutralization  scenarios,  being  coupled  with  empirical  models  which  provide  the 
initial  post-rupture  dispersed  phase  conditions.  Emphasis  is  on  biological  agent  scenarios  and 
detailed  modeling  of  the  spore  is  required  to  obtain  an  accurate  portrayal  of  how  thermal  energy  on 
the  surface  is  transferred  internally  (this  goes  well  beyond  the  need  to  include  T(r)  effects  in  liquid 
droplets).  Other  work  with  CRAFT  in  this  arena  entails  investigation  of  post-hit  cruise  missile 
scenarios  using  high  energy  mechanisms  to  neutralize  agents. 

The  SDROP  work  was  performed  as  a  side  activity  to  provide  guidance  to  Army  and  Navy 
groups  relevant  to  high  altitude  particle  fate  simulations.  As  described,  higher-order  modeling 
details  such  as  droplet  deformation  were  shown  to  have  a  first-order  effect  on  particle  fate  and 
cannot  be  neglected. 
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APPENDIX  A 


Atmospheric  Properties 


Using  the  results  of  (Ref.  24),  the  temperature  of  the  atmosphere  is  eomputed  as: 


M  n 

jrW 

Mq  0 


(A.l) 


In  the  current  implementation  of  SDROP,  n  =  8,  and  the  coefficients  are  given  by: 


Oo 

2.824791081e2 

a5 

6.050186406e-6 

-5.240572992 

“6 

-3.550162735e-8 

02 

-1.266010595e-l 

1.014102927e-10 

o3 

1.873293836e-2 

as 

-1.124449619e-13 

o4 

-5.104746533e-4 

For  Z  <  90  [km} ,  the  molecular  weight  is  constant: 


M=M0=28.9644 


(A.2) 


whereas,  above  90  [km],  a  5th-order  polynomial  was  fit  to  the  data  of  (Ref.  24)  to  obtain: 

The  pressure  and  density  are  computed  from: 

-£"  =  exp  (A.3) 

*0  K  lm 

-^•  =  ^-exp  (~M±\lJLdZ)  (A.4) 

Po  lM  K  IM 

where  the  integral  is  evaluated  as: 

\l~~dZ  =  b0x  [In (z+b3  )]o  -  etc.  (A.5) 

*  M 


with  the  coefficients: 
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*0 

-1.0902039e-7 

b\2 

-1.60023  le-  4 

bi 

6356.77 

b\z 

-189.5201 

b2 

1.787026^-4 

bu 

9665.295 

21.680485 

bl5 

1.163707k -3 
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-1.3949832* -5 

^16 
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-284.01768 

b\7 
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1.3327563* -4 

b\z 

-5.562892* -5 

-29.89506 

bl9 

-420.11368 

b& 
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b2o 

45,675.466 

b9 

8.3 168074* -4 

b2l 

1.284404e-4 

b\o 

0.037777365 

b22 

0.025387008 

b\\ 

0.5646783 

bn 

-5.3327146 

The  atmospheric  viscosity  is  computed  as: 


bT r 
T+S 


(A.6) 


where  b  =  1.458e-6  [kg/smK05]  and  S=1 10.4  [K]  is  Sutherland’s  constant.  The  above  expression, 
which  is  derived  from  kinetic  theory  and  fit  to  experiments,  “...fails  for  conditions  of  very  high  and 
very  low  temperatures  and  under  conditions  occurring  at  great  altitudes.  (Consequently  tabular 
entries  for  coefficient  of  viscosity  are  terminated  at  Z  =  90  [km].)” 
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APPENDIX  B 


Liquid  Phase  Thermodynamic  Properties 
B-l.  Agent  and  Agent  Simulants 

The  thermodynamic  properties  of  chemical  agents  and  agent  simulants  are  computed  using 
the  approximate  formulae  from  the  VLSTRACK  code,  which  are  given  in  cgs  units.  Accordingly, 
the  liquid  density  is: 

Pp  =  Pa~  Pbx  [#]  (B.l) 

where  TpC  is  the  drop  temperature  in  Celsius.  The  specific  enthalpy  of  vaporization  is: 

T 

f 

^  i  +c 

A Hy  =  2.303 RB — - 4.184  x  107  [erg/  g]  (B.2) 

Mp 

The  hquid  viscosity  is: 

(HA  + — — — ) 

Pp  =  10"2  x  pp  x  10  Tp-c+Mc  [ g/cms ]  (B.3) 

The  diffusivity  of  simulant  vapor  into  air  is: 


4 


t 


D  =4.3x10 


-2 


[v°3333+( 


0-3333  j2 

P? 


[cm2  /  s] 


(B.4) 


The  surface  tension  is: 

Gp  =  oA-  GBTpC  [dyn/ cm]  (B.5) 

* 

Finally,  the  vapor  pressure  of  the  simulant  is  obtained  by  applying  the  perfect  gas  law  to  the 
approximate  formulae  for  the  volatility,  which  is  defined  as  the  density  of  vapor  at  saturation 
conditions. 

Pp>v  =  1 .6  x  10"510  TpC+c  R  (B.6) 

For  VX,  the  constants  used  in  the  above  expressions  are: 


-79- 


Pa 

1.029 

Pc 

80 

Pb 

8.3 

VA 

29.9 

A 

7.281 

<*A 

34.6631 

B 

2072.1 

0.1326 

C 

172.5 

Pa 

29 

Pa 

-0.80572 

_ Bb 

189.38 

R 

1.9872 

Remaining  properties  for  VX  that  were  obtained  from  DROPFATE  and  VLSTRACK  people  are: 


Pc[g/cc] 

0.3177 

K[W/mK] 

0.115 

Pcfatm] 

19.0 

CP[cal/gC] 

0.447  @  35  [C] 

Tr[K] 

746.43 

Tb[c] 

267.4 

Vr[cc/mol] 

841.66 

Tf[C] 

-50.0 

B-2.  Other  Liquids 

For  liquids  other  than  agent  and  agent  simulants,  the  thermodynamic  properties  were 
computed  using  approximate  formulae  and  rules  from  the  reference  work  of  Reid  [25].  The  vapor 
pressure  is  computed  using  the  Wagner  equation: 


p„,  ■  expt-^ (1  - ^ + W>>(1  - r')''5 ±VMzlf± YM - Lt )Pr  ll0sPa]  (B.7) 

The  specific  enthalpy  of  vaporization  is  computed  with  a  combination  of  Vetere  method  and 
Fish  and  Lielmezs  method.  The  former  is  used  to  compute  the  enthalpy  of  vaporization  at  the 
normal  boiling  point: 


Wvb=R7 'cTbr 


0.4343Itt(Pc)  -  0.6943 1  +  0.895847;, 
0.37691  -  0.373067],,.  +  0.15075P^‘7r2 


l J/kg ] 


(B.8) 


'C  *br 


and  the  latter  provides  the  necessary  temperature  correction: 


A#v=A Hvb 


T,  X  +  X? 
Tbr  l  +  Xp 


[J/kg], 


X  =  IkL  \z!l 


Tr  1  ~Tbr 


(B.9) 


For  inorganic  and  organic  liquids,  q  =  0.35298  and  p  =  0.13856. 

The  Cp  of  the  liquid  is  obtained  as  a  correction  to  the  Cp  of  the  pure  gas  in  the  ideal-gas 
state  using  the  corresponding  state  method  of  Rowlinson  and  Bondi: 
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c,.,+cj. 


=  1.45  +  0.45(1  -  Tr)~x  +0.25® [17.1 1  +  25.2(1  -  J,)3  271  +1.742(1  -  Tr)~x ]  (B.10) 


The  thermal  conductivity  of  organic  liquids  may  be  computed  using  the  boiling-point 
method  of  Sato  to  compute  the  conductivity  at  the  normal  boiling  point  and  the  Riedel  equation  as  a 
temperature  correction: 


2 


1.11  [3  +  20(l-rr)3] 

1  2  * 


M2  [3  +  20(l-7],r)3] 


[' W/mK ] 


(B.ll) 


Note  that  this  and  other  approximations  are  valid  only  in  the  temperature  region  below  the  normal 
boiling  point,  where  the  thermal  conductivity  does  not  vary  greatly.  From  Reid  [25],  use  of  this 
relationship  for  Tr  >  0.65  is  not  recommended. 


APPENDIX  B 


-81  - 


